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Introduction 



In which we introduce some of the basic concepts of modern density functional theory, 
including the Kohn-Sham description of a system, and give a simple but powerful example 
of DFT at work. 

1.1 Importance 

1 Density functional theory (DFT) has long been the mainstay of electronic structure calcula- 
tions in solid-state physics. In the 19990's it became very popular in quantum chemistry. This 
is because approximate functionals were shown to provide a useful balance between accuracy 
and computational cost. This allowed much larger systems to be treated than by traditional 
ab initio methods, while retaining much of their accuracy. Nowadays, traditional wavefunc- 
tion methods, either variational or perturbative, can be applied to find highly accurate results 
on smaller systems, providing benchmarks for developing density functionals, which can then 
be applied to much larger systems 

But DFT is not just another way of solving the Schrodinger equation. Nor is it simply a 
method of parametrizing empirical results. Density functional theory is a completely different, 
formally rigorous, way of approaching any interacting problem, by mapping it exactly to a 
much easier-to-solve non-interacting problem. Its methodology is applied in a large variety of 
fields to many different problems, with the ground-state electronic structure problem simply 
being the most common. 

The aim of this book is to provide a relatively gentle, but nonetheless rigorous, introduction 
to this subject. The technical level is no higher than any graduate quantum course, or many 
advanced undergraduate courses, but leaps at the conceptual level are required. These are 
almost as large as those in going from classical to quantum mechanics. In some sense, they 
are more difficult, as these leaps are usually made when we are more advanced (i.e., set) in 
our thinking. 

Students from all areas of modern computational science (chemistry, physics, materials 
science, biochemistry, geophysics, etc.) are invited to work through the material. The only 

1 ©2000 by Kieron Burke. All rights reserved. 
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necessary requirements are a good background in elementary quantum mechanics, no fear 
of calculus of more than one variable, and a desire to learn. The student should end up 
knowing what the Kohn-Sham equations are, what functionals are, how much (or little) is 
known of their exact properties, how they can be approximated, and how insight into all these 
things produces understanding of the errors in electronic structure calculations. Preliminary 
forms of my lecture notes have been used throughout the world during 5 years its taken to 
complete this book. The hope is that these notes will be used by students worldwide to gain 
a better understanding of this fundamental theory. I ask only that you send me an email (to 
kieron@rutchem.rutgers.edu) if you use this material. In return, I am happy to grade problems 
and answer questions for all who are interested. These notes are copyright of Kieron Burke. 
No reproduction for purposes of sale is allowed. 

The book is laid out in the form of an undergraduate text, and requires the working of 
many exercises (although many of the answers are given). The idea is that the book and 
exercises should be easy reading, but leave the reader with very clear concepts of modern 
density functional theory. Throughout the text, there are exercises that must be performed 
to get full value from the book. Also, at the end of each chapter, there are questions aimed 
at making you think about the material. These should be thought about and answere as you 
go along, but answers don't need to be written out as explicitly as for the problems. 

1.2 What is a Kohn-Sham calculation? 

To give an idea of what DFT is all about, and why it is so useful, we start with a very simple 
example, the hydrogen molecule, Hi- Throughout this book, we make the Born-Oppenheimer 
approximation, in which we treat the heavy nuclei as fixed points, and we want only to solve 
the ground-state quantum mechanical problem for the electrons. 

In regular quantum mechanics, we must solve the interacting Schrodinger equation: 

-\ £ V 2 + ^^ + £ «er t (r i )U(ri,r 2 ) = S*(r 1 ,r 2 ), (1.1) 
2j=i,2 |ri-r 2 | i= L2 J 

where the index % runs over the two electrons, and the external potential, the potential 
experienced by the electrons due to the nuclei, is 

v ext (r) = -Z/r-Z/\r-Rz\, (1.2) 

where Z = 1 is the charge on each nucleus, z is a unit vector along the bond axis, and R is 
a chosen internuclear separation. Except where noted, we use atomic units throughout this 
text, so that 

e 2 = h = m, = 1, (1.3) 

where e is the electronic charge, % is Planck's constant, and m is the electronic mass. As 
a consequence, all energies are in Hartrees (1 H = 27.2114 eV= 627.5 kcal/mol) and all 
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distances are given in Bohr radii (a = 0.529A). In this example, the electrons are in a 
spin singlet, so that their spatial wavefunction $(ri,r 2 ) is symmetric under interchange of 
ri and r 2 . Solution of Eq. (1.1) is complicated by the electrostatic repulsion between the 
particles, which we denote as V ee . It couples the two coordinates together, making Eq. (1.1) 
a complicated partial differential equation in 6 coordinates, and its exact solution can be quite 
demanding. In Fig. 1.1 we plot the results of such a calculation for the total energy of the 
molecule, E + l/R, the second term being the Coulomb repulsion of the nuclei. The position 
of the minimum is the equilibrium bond length, while the depth of the minimum, minus the 
zero point vibrational energy, is the bond energy. More generally, the global energy minimum 
determines all the geometry of a molecule, or the lattice structure of a solid, as well as all the 
vibrations and rotations. But for larger systems with N electrons, the wavefunction depends 
on all 3N coordinates of those electrons. 




R(A) 

Figure 1.1: The total energy of the H2 molecule as a function of internuclear separation. 

We note at this point that, with an exact ground-state wavefunction, it is easy to calculate 
the probability density of the system: 

n(r) =2/dV |*(r,r')| 2 . (1.4) 

The probability density tells you that the probability of finding an electron in d 3 r around r is 
n{v)d?r. For our H 2 molecule at equilibrium, this would look like the familiar two decaying 
exponentials centered over the nuceli, with an enhancement in between, where the chemical 
bond has formed. 

Next, imagine a system of two non-interacting electrons in some potential, v s (r), chosen 
somehow to mimic the true electronic system. Because the electrons are non-interacting, 
their coordinates decouple, and their wavefunction is a simple product of one-electron wave- 
functions, called orbitals, satisfying: 

{-\v 2 + v s (r)^ i (v) = e i ^l (1-5) 
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where $(r 1 ,r 2 ) = 0o(r 1 )0 o (r 2 ). This is a much simpler set of equations to solve, since it 
only has 3 coordinates. Even with many electrons, say N, one would still need to solve only 
a 3-D equation, and then occupy the first N/2 levels, as opposed to solving a 3iV-coordinate 
Schrodinger equation. If we can get our non-interacting system to accurately 'mimic' the 
true system, then we will have a computationally much more tractable problem to solve. 

How do we get this mimicking? Traditionally, if we think of approximating the true wave- 
function by a non-interacting product of orbitals, and then minimize the energy, we find the 
Hartree-Fock equations, which yield an effective potential: 2 

« s HF (r)=« ext (r) + i/dV^. (1.6) 

The correction to the external potential mimics the effect of the second electron, in particular 
screening the nuclei. For example, at large distances from the molecule, this potential decays 
as — 1/r, reflecting an effective charge of Z — 1. Note that insertion of this potential into 
Eq. (1.5) now yields a potential that depends on the electronic density, which in turn is 
calculated from the solution to the equation. This is termed therefore a self-consistent set 
of equations. An initial guess might be made for the potential, the eigenvalue problem is 
then solved, the density calculated, and a new potential found. These steps are repeated 
until there is no change in the output from one cycle to the next - self-consistency has been 
reached. Such a set of equations are often called self-consistent field (SCF) equations. In 
Fig. 1.1, we plot the Hartree-Fock result and find that, although its minimum position is very 
accurate, it underbinds the molecule significantly. This has been a well-known deficiency of 
this method, and traditional methods attempt to improve the wavefunction to get a better 
energy. The missing piece of energy is called the correlation energy. 

In a Kohn-Sham calculation, the basic steps are very much the same, but the logic is 
entirely different. Imagine a pair of non-interacting electrons which have precisely the same 
density n(r) as the physical system. This is the Kohn-Sham system, and using density 
functional methods, one can derive its potential v s (r) if one knows how the total energy 
E depends on the density. A single simple approximation for the unknown dependence of 
the energy on the density can be applied to all electronic systems, and predicts both the 
energy and the self-consistent potential for the fictitious non-interacting electrons. In this 
view, the Kohn-Sham wavefunction of orbitals is not considered an approximation to the 
exact wavefunction. Rather it is a precisely-defined property of any electronic system, which 
is determined uniquely by the density. To emphasize this point, consider our H 2 example in 
the united atom limit, i.e., He. In Fig. 1.2, a highly accurate many-body wavefunction for 
the He atom was calculated, and the density extracted. In the bottom of the figure, we plot 
both the physical external potential, —2/r, and the exact Kohn-Sham potential. 3 Two non- 
interacting electrons sitting in this potential have precisely the same density as the interacting 

2 For the well-informed, we note that the correction to the external potential consists of the Hartree potential, which is double that shown above, 
less the exchange potential, which in this case cancels exactly half the Hartree. 

3 It is a simple exercise to extract the exact Kohn-Sham potential from the exact density, as in section ??. 
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Figure 1.2: The external and Kohn-Sham potentials for the He atom, thanks to Cyrus' Umrigar's very accurate density. 

electrons. If we can figure out some way to approximate this potential accurately, we have 
a much less demanding set of equations to solve than those of the true system. Thus we 
are always trying to improve a non-interacting calculation of a non-interacting wavefunction, 
rather than that of the full physical system. In Fig. 1.1 there are also plotted the local 
density approximation (LDA) and generalized gradient approximation (GGA) curves. LDA 
is the simplest possible density functional approximation, and it already greatly improves on 
HF, although it typically overbinds by about 1/20 of a Hartree (or 1 eV or 30 kcal/mol), 
which is too inaccurate for most quantum chemical purposes, but sufficiently reliable for many 
solid-state calculations. More sophisticated GGA's (and hybrids) reduce the typical error in 
LDA by about a factor of 5 (or more), making DFT a very useful tool in quantum chemistry. 
In section 1.5, we show how density functionals work with a simple example from elementary 
quantum mechanics. 

1.3 Reinterpreting molecular orbitals 

The process of bonding between molecules is shown in introductory chemistry textbooks as 
linear combinations of atomic orbitals forming molecular orbitals of lower energy, as in Fig. 
1.3. 

But later, in studying computational chemistry, we discover this is only the Hartree-Fock 
picture, which, as stated above, is rarely accurate enough for quantum chemical calculations. 
In this picture, we need a more accurate wavefunction, but then lose this simple picture of 
chemical bonding. This is a paradox, as chemical reactivity is usually thought of in terms of 
frontier orbitals. 

In the Kohn-Sham approach, the orbitals are exact and unique, i.e., there exists (at most) 
one external potential that, when doubly occupied by two non-interacting electrons, yields 
the exact density of the H 2 molecule. So in this view, molecular orbital pictures retain their 
significance, if they are the exact Kohn-Sham orbitals, rather than those of Hartree-Fock. 
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Figure 1.3: Orbital diagram of H2 bond formation. 

And a highly accurate approximate density functional calculation produces the full electronic 
energy from these orbitals, resolving the paradox. 

1.4 Sampler of applications 

In this section, we highlight a few recent applications of modern DFT, to give the reader a 
feeling for what kinds of systems can be tackled. 

Example from biochemistry, using ONIOM. 

Example from solid-state, eg ferromagnetism. 

Example using CP molecular dynamics, eg phase-diagram of C. 

Example from photochemistry, using TDDFT. 

1.5 Particle in a box 

In this section, we take the simplest example from elementary quantum mechanics, and apply 
density functional theory to it. This provides the underlying concepts behind what is going 
on in the more previous sections of this chapter. Much of the rest of the book is spent 
connecting these two. 

In general, we write our Hamiltonian as 

H = f + V (1.7) 

where T denotes the operator for the kinetic energy and V the potential energy. We begin 
with the simplest possible case. The Hamiltonian for a 1-D 1-electron system can be written 
as 

1 d 2 ._. . . 
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The time-independent Schrodinger equation has solutions: 

Hfa = etfi, i = 0,l,... (1.9) 

Thus e denotes the ground state energy and <p the ground state wave function. Because 
the operator H is hermitian the eigenstates can be chosen orthonormal: 

dxfa(x)fa(x) = 5 i j (1.10) 
Thus each eigenfunction is normalized. The electron probability density is just n[x) = 

l<M*)l 2 . 

Exercise 1 Particle in a one-dimensional box 

Consider the elementary example of a particle in a box, i.e., V(x) = oc everywhere, except 
< x < L, where V = 0. This problem is given in just about all elementary textbooks on 
quantum mechanics. Show that the solution consists of trigonometric functions (standing 
waves), and making them vanish at the boundaries quantizes the energy, i.e., 

& = ^sin(fcix), i = l,2,... (1.11) 

where fc t = m/L, and the energies are e, = kj/2. Check the orthonormality condition, Eq. 
(1.10). 

Next we consider the following approximate density functional for the kinetic energy of 
non-interacting electrons in one-dimension. Don't worry if you never heard of a functional 
before. In fact, chapter ?? discusses functionals in some detail. For now, we need only the 
fact that a functional is a rule which maps a function onto a number. We use square brackets 
[..] to indicate a functional dependence. The functional we write down will be the local 
approximation to the kinetic energy of non-interacting spinless fermions in one dimension. 
We do not derive it here, but present it for calculational use, and derive it later in the chapter: 

T loc [n] = 1.645 l°° dx n 3 (x). (1.12) 

A local functional is one which is a simple integral over a function of its argument. Now, 
since for the particle in a box, all its energy is kinetic, we simply estimate the energy using 
T s ! oc [rt]. Since the electron density is the square of the orbital: 

2 

ni (a;) = -sin 2 (fc 1 x), (1-13) 

Li 

we find T s loc = 4.11/L 2 , a 17% underestimate of the true value, tt 2 /2L 2 = 4.93/L 2 . What 
is so great about that? There is not much more work in finding the exact energy (just take 
the second derivative of fa, and divide by fa) as there is in evaluating the integral. 

But watch what happens if there are two particles, fermions of the same spin (the general 
case for which T s loc was designed). We find now E = e, + e 2 = 5n 2 /2L 2 = 24.7/L 2 . If we 
evaluate the approximate kinetic energy again, using n(x) = n\(x) + j j sm 2 (k2x), we find 
r s loc = 21.8/L 2 , an 11% underestimate. 
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Exercise 2 Three particles in a box 

Calculate the exact energy of three identical fermions in a box of length L. Plot the total 
density for one, two, and three particles. Compare the exact energy with that from the local 
approximation. What happens as the number of electrons grows? Answer: E = 7ir 2 /L 2 , 
E [oc = 63/L 2 , a 9% underestimate. 

These results are gotten by the evaluation of a simple integral over the density, whereas 
the exact numbers require solving for all the eigenvalues of the box, up to N, the number of 
electrons in it. Thus approximate density functionals produce reliable but inexact results at 
a fraction of the usual cost. 

To see how remarkable this is, consider another paradigm of textbook quantum mechanics, 
namely the one-dimensional harmonic oscillator. In the previous example, the particles were 
sine waves; in this one, they are Gaussians, i.e., of a completely different shape. Applying 
the same approximate functional, we find merely a 20% overestimate for one particle: 

Exercise 3 Kinetic energy of 1-d harmonic oscillator 

What is the ground-state energy of a 1-d harmonic oscillator? What is its kinetic energy? 
Calculate its density, and from it extract the local density approximation to its kinetic energy. 
Repeat for two same-spin electrons occupying the lowest two levels of the well. 

Congratulations. You have performed your first elementary density functional calculations, 
and are well on the way to becoming an expert! 
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Figure 1.4: Density of 10 particles in a box. 

In fact, there's a simple way to see how this magic has been performed. Imagine the box 
with a large number of particles. In Fig. 1.4, we have plotted the case for N = 10. As the 
number of particles grows, the density becomes more and more uniform (independent of the 
shape of the box!). If we assume it actually becomes uniform as N — > oo, with errors of 
order 1/N or smaller, this gives us the constant in the local approximation. To see this, note 
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that the exact energy is 

E=^ 2 T,f = ^N(N + l)(2N + l)/6 (1.14) 

whose leading contribution is simply 7r 2 iV 3 /6. Since the density has become uniform (to this 
order), the local approximation becomes exact, yielding the constant. 

Later in the book, we will see how, if we had an extremely accurate non-interacting kinetic 
energy density functional in three dimensions, we could revolutionize electronic structure 
calculations, by making them very fast. This is because we would then have a method for 
calculating the density and ground-state energy of systems which involved solving a single 
self-consistent integro-differential equation for any electronic system. This would be a true 
density functional calculation. Most present calculations do not do this. By solving an 
effective single-particle equation for the orbitals (the Kohn-Sham equation), they find T s 
exactly, but at large computational cost (for large systems). So now you even know an 
unsolved problem in density functional theory. 

Now answer the following simple questions as well as you can, justifying your answers in 
each case. 

1.6 Questions 

We begin with conceptual questions, that are designed to test your understanding of the 
material presented in the chapter. 

1. Why is a Kohn-Sham calculation much faster than a traditional wavefunction calculation? 

2. If you evaluate the kinetic energy of a Kohn-Sham system, is it equal to the physical 
kinetic energy? 

3. Repeat above question for (1/r), which can be measured in scattering experiments. 

4. Why is the density cubed in the local approximation for T s ? (see section ?? for the 
answer). 

We continue with more creative questions, designed to make you think about what is in 
the chapter. 

5. Suppose you'd been told that T s loc is proportional to / dxn 3 (x), but were not told the 
constant of proportionality. If you can do any calculation you like, what procedure might 
you use to determine the constant in Tj| oc ? (see section 11.5 for the answer). 

6. In what way will T^° c change if spin is included, e.g., for two electrons of opposite spin 
in a box? (see section 9.2 for the answer). 

7. If we add an infinitesimal to n(x) at a point, i.e., e5(x) as e — > 0, how does T s loc change? 
(see section 2.2 for the answer). 
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8. The simplest density functional approximation is a local one. What form might a cor- 
rection to Tl° c take? (see section 17.2 for the answer). 

9. In the chapter and exercises, T s loc does very well. But there are cases where it fails quite 
badly. Can you find one, and say why? 

10. Consider a single electron in an excited state, with density n*(x). Will T s loc [n*] provide 
an accurate estimate of its kinetic energy? 
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Extra exercises 

1. For a J-function potential, v(x) = —ZS(x), the ground-state kinetic energy is Z 2 /2 and 
the density is Z exp(—2Zx). Before calculating the local approximation, first guess if it 
will be an overestimate or an underestimate, and by how much. 

2. Another famous example is v(x) = — a 2 sech 2 (ax), with ground-state energy — a 2 /2 and 
density asech 2 (ox)/2. Repeat above problem. 



Chapter 2 
Functionals 



1 In this chapter, we introduce in a more systematic fashion what exactly a functional 
is, and how to perform elementary operations on functionals. This mathematics will 
be needed when we describe the quantum mechanics of interacting electrons in terms of 
density functional theory. 

2.1 What is a functional? 

A function maps one number to another. A functional assigns a number to a function. For 
example, consider all functions r(9), < 9 < 2ir, which are periodic, i.e., r(9 + 2ir) = r(9). 
Such functions describe shapes in two-dimensions of curves which do not "double-back" on 
themselves, such as in Fig. 2.1. For every such curve, we can define the perimeter P as the 



Figure 2.1: A 2D curve which is generated by a function r — r(0) 

length of the curve, and the area A as the area enclosed by it. These are functionals of r (6), 
in the sense that, for a given curve, such as the ellipse 

r{9) = l/yW^) + 4cos 2 (6>) (2.1) 

there is a single well-defined value of P and of A. We write P[r] and A[r] to indicate this 
functional dependence. Note that, even if we don't know the relation explicitly, we do know 
it exists: Every bounded curve has a perimeter and an area. 

For this simple example, we can use elementary trigonometry to deduce explicit formulas 
for these functionals. Consider the contribution from an infinitesimal change in angle d6, and 

1 ©2000 by Kieron Burke. All rights reserved. 
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integrate over the entire range of angles. The contribution to the area is that of a triangle 
of base rcos(dff) « r and height r(0)sm.(dO) « rd6, yielding 

1 r 1 r2rr 

A[r] = -f rd9r = -J q d9 r 2 {9) (2.2) 

For an abitrary r(9) the integral above can be evaluated to give the area that is enclosed by 
the curve. Thus this functional maps a real function of one argument to a number, in this 
case the area enclosed by the curve. Similarly, the infinitesimal change in the perimeter is 
just the line segment in polar coordinates: 

SP 2 = dr 2 + r 2 {9)d9. (2.3) 

Now, since r = r(ff), we can write dr = d9 (dr/dff), yielding 

P[r] = J*"d6 <Jr 2 {9) + (dr/d9) 2 . (2.4) 

In both cases, once we know r(9), we can calculate the functional. 

The area is a local functional of r(0), since it can be written in the form: 

A[r] = f^d9f(r(9)), (2.5) 

where f(r) is a function of r. For the area, / = r 2 /2. These are called local because, 
inside the integral, one needs only to know the function right at a single point to evaluate 
the contribution to the functional from that point. On the other hand, P is a semi-local 
functional of the radius, as it depends not only on r, but dr/d9. 

We have already come across a few examples of density functionals. In the opening 
chapter, the local approximation to the kinetic energy is a local functional of the density, 
with f(n) = 1.645n 3 . On the other hand, for any one-electron system, the exact kinetic 
energy is a semi-local functional, called the von Weisacker functional: 

rVWM = j£>£ (2-6) 
where n'(x) = dn/dx. Later we will see some examples of fully non-local functionals. 

2.2 Functional derivatives 

When we show that the ground-state energy of a quantum mechanical system is a functional 
of the density, we will then want to minimize that energy to find the true ground-state density. 
To do this, we must learn how to differentiate functionals, in much the same way as we learn 
how to differentiate regular functions in elementary calculus. 

To begin with, we must define a functional derivative. Imagine making a tiny increase in 
a function, localized to one point, and asking how the value of a functional has changed due 
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to this increase, i.e., we add an infinitesimal change 5r(9) = e5(Q — 9 ). 2 How does A[r\ 
change? 

A[r + 5r] - A[r] = * J** d9 {(r + e5(6 - 9 )f - r 2 } = d9 re5{9 - 9 ) = e r(9 ) (2.7) 

The functional derivative, denoted SA/5r(9) is just the change in A divided by e, or just r(9) 
in this case. Since this is linear in the change in r, the general definition, for any infinitesimal 
change in r, is just 

A[r + 8r]-A[r] = J Q d9 ^Sr(8), (2.8) 

where the functional derivative 5A/Sr(8) is that function of 9 which makes this formula 
exact for any small change in r{9). Just as the usual derivative df /dx of a function f(x) 
tells you how much / changes when x changes by a small amount, i.e., f(x + dx) — f(x) = 
(df/dx) dx + 0(dx 2 ), so does the functional derivative, a function, tell you how much a 
functional changes when its arguments changes by a small "amount" (in this case, a small 
function). Thus 

VA[rm = W) =r{0) ' (2 ' 9) 
where the notation va[t](9) emphasizes that the functional derivative of A[r] is a ^-dependent 
functional. 

Exercise 4 Derivative of a local functional 

Show that, for any local functional A[n] = J dxa(n(x j), the functional derivative is S A/8n(x) - 
a'(n(x)), where a'(n) = da/dn. 

In general, the way to find a functional derivative is to evaluate the expression A[r+(5r] — A[r] 
to leading order in 5r, and the resulting integral must be cast in the form of a function times 
8r. Often it is necessary to do an integration by parts to identify the functional derivative. 
For the perimeter, such complications arise. We have 

P[r + Sr] = d8 J(r + Sr) 2 + {? + 5r' f (2.10) 

where r' = dr/dO and all functions are assumed to have argument 9. To first order in 5r, we 
find 

P[r + Sr] = fd9^^(l + r % + + r ^ 

= p M + f de Sr + 5r ') (2 ' n) 



2 For purists, we take a Gaussian of tiny but finite width, so that as e — > 0, the maximum change becomes arbitrarily small. Then we take the 
width to zero! 
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The first integral is in just the right form for identifying the contribution to the functional 
derivative, but the second is not. But we can perform an integration by parts to find 

5P r d ( r' \ 

frrJO) = P[rm = vV 2 + r' 2 " M [Vr* + r' 2 j ^ 

Note that the end-point term of the integration-by-parts vanishes, because r(9) is periodic. 
A little calculus and algebra finally yields 

„3 I Orr' 2 — r 2 r") 
Exercise 5 Derivative of a semilocal functional 

Show that the functional derivative of a semilocal functional B[n] = J^dx b{n,dn/dx), 
assuming n and its derivatives vanish rapidly as \x\ — > oc, is 

r t , x dh d I db\ ,„ , 

MnlW^-T-b ( 2 - 14 ) 



dn dx \dn 
where n'(x) = dn/dx. Use your answer to show 

fiJiVW „!l 



n n 

+ (2-15) 



Sn(x) An 8n 2 
The density functional we use in practice are three-dimensional: 

Exercise 6 Hartree potential 

The Hartree (or classical electrostatic) energy of a charge distribution interacting with itself 
via Coulomb 's law is given by 



U [n] = y*rj*r>'^. (2.16) 
Show that its functional derivative, the Hartree potential, is 

v "["]W = / rfV (S- ( 2 - 17 ) 



2.3 Euler-Lagrange equations 

The last piece of functional technology we need for now is how to solve a constrained opti- 
mization problem, i.e., how to maximize (or minimize) one functional, subject to a constraint 
imposed by another. 

For example, we might want to know what is the maximum area we can enclose inside a 
loop of string of fixed length /. Thus we need to maximize A[r], subject to the constraint 
that P[r] = I. There is a well-known method for doing such problems, called the method 
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of Lagrange multipliers. Following the example in the Appendix (A.l), construct a new 
functional 

B[r]=A[r]- l iP[r], (2.18) 

where fj, is at this point an unknown constant. Then extremize the new functional B[r], by 
setting its functional derivative to zero: 

5B 5 A SP r i + r' 2 (2r-r") „ 

* = r{0) ~ il +>)»/» = °' ( 2 - 19 ) 

A simple solution to this equation is r(6) = /i, a constant. This tells us that the optimum 
shape is a circle, and the radius of that circle can be found by inserting the solution in the 
constraint, P[r = p] = 2np, = I. The largest area enclosable by a piece of string of length / 

is A[r = n\ = Z 2 /(4tt). 

Exercise 7 Change in a functional 

Suppose you know that a functional E[n] has functional derivative v(r) = —l/r for the 
density n(r) = Z exp(—2Zr) , where Z = 1. Estimate the change in E when Z becomes 
1.1. 

Exercise 8 Second functional derivative 

Find the second functional derivative of (a) a local functional and (b) T s vw . 



2.4 Questions 

1. Compare the functional derivative of T/ W [n] with T g loc [n] for some sample one-electron 
problem. Comment. 

2. If someone just tells you a number for any density you give them, e.g., the someone might 
be Mother Nature, and the number might be the total energy measured by experiment, 
devise a method for deducing if Mother Nature's functional is local or not. 

3. Is there a simple relationship between T s and J dx n(x)ST s /5n(x)7 First consider the 
local approximation, then the Von Weisacker. Comment on your result. 

4. For fixed particle number, is there any indeterminancy in the functional derivative of a 
density functional? 
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Extra exercises 

1. Consider a square of side 1. Check that the circle of the same perimeter has a larger 
area. 

2. For the unit square centered on the origin, check that Eqs. (2.2) and (2.4). are correct. 
Hint: You need only do the first octant, x > y > 0. 



Chapter 3 



One electron 



1 In this chapter, we review the traditional wavefunction picture of Schrddinger for one 
particle, but introducing our own specific notation. 

3.1 Variational principle 

We start off with the simplest possible case, one electron in one dimension. Recall, from 
basic quantum mechanics, the Rayleigh-Ritz variational principle: 

E = mm(<j)\H\0), j°° dx U(x)\ 2 = 1. (3.1) 

<p J— CO 

We will use this basic, extremely powerful principle throughout this book. This principle says 
that we can use any normalized wavefunction to calculate the expectation value of the energy 
for our problem, and we are guaranteed to get an energy above the true ground-state energy. 
For example, we can use the same wavefunction for every 1-d one-electron problem, and get 
an upper bound on the ground-state energy. 

Having learnt about functionals in Chapter ??, we may write this principle in the following 
useful functional form. For example, the kinetic energy of the particle can be considered as 
a functional of the wavefunction: 

T[4>] = (-\^ fa) = \ £>K(*)| 2 (3.2) 

where the second form is gotten by integration by parts, and the prime denotes a spatial 
derivative. (This second form is much handier for many calculations, as you only need take 
one derivative.) Thus, for any given normalized wavefunction, there is a single number T, 
which can be calculated from it, via Eq. (3.2). Similarly, the potential energy is a very simple 
functional: 

V[</>] = /_^dxl/(z) !<?(*) | 2 . (3.3) 

The kinetic energy of a given wavefunction is always the same, no matter what the prob- 
lem we are applying that wavefunction to. The kinetic energy is a universal functional, 

1 ©2000 by Kieron Burke. All rights reserved. 
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independent of the particular problem, i.e., we apply the same operation on a given trial 
wavefunction, no matter what our physical problem is. But the potential energy functional 
differs in each problem; it is not universal. This is analogous to trying to find the minimum 
over x of f(x,y) = x 2 /2 — xy. The first term is independent of y, so we could make a list 
of it value for every x, and then use that same list to find the minimum for any y, without 
having to evaluate it again for each value of y. 

Now our variational principle can be written by writing 

E = mm{T[d>\ + V[0\}, f° dx |0(.x)| 2 = 1 (3.4) 

We need simply evaluate the energy of all possible normalized wavefunctions, and choose the 
lowest one. 

3.2 Trial wavefunctions 

Next, we put trial wavefunctions into the variational principle to find upper bounds on the 
ground-state energy of a system. 

Exercise 9 Gaussian trial wavefunction 

An interesting potential in one-dimension is just V(x) = —5{x), where 6 is the Dirac delta 
function. Using a Gaussian wavefunction, 4>c,{x) = n~ l l A exp(— x 2 /2), as a trial wavefunc- 
tion, make an estimate for the energy, and a rigorous statement about the exact energy. 

We can improve on our previous answers in two distinct ways. The first method is to 
include an adjustable parameter, e.g., the length scale, into our trial wavefunction. 

Exercise 10 Adjustable trial wavefunction 

Repeat Ex. (9) with a trial wavefunction <pc(ax) (don't forget to renormalize). Find (T)(ct) 
and ( V )(ct) and plot their sum as a function of a. Which value of a is the best, and what 
is your estimate of the ground-state energy? Compare with previous calculation. 

Note that this has led to a non-linear optimization problem in a. Varying exponents in 
trial wavefunctions leads to difficult optimization problems. 

To see how well we did for this problem, the next exercise yields the exact answer. 

Exercise 11 1-d Hydrogen atom 

Use (t> E (ax), where 4>e{x) = exp(— \x\) as a trial wavefunction for the problem in Ex. (9). 
Find the lowest energy. This is the exact answer. Calculate the errors made in the previous 
exercises, and comment. 

An alternative, simpler approach to improving our trial wavefunction is to make a linear 
combination with other trial wavefunctions, but not vary the length scale. Thus we write 



<t>trial{x) = C (t> {x) + Ci</i>i(z), 



(3.5) 
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and minimize the expectation value of H, subject to the constraint that <j> be normalized. 
This leads to a set of linear equations 

||H-BS||=0, (3.6) 

where 

/oo 
dx <j>*i{x)H4>j{x) (3.7) 

is the Hamiltonian matrix, and 

/oo 
^dx 4>*{x)4>j{x) (3.8) 

is the overlap matrix. In general, Eq. (3.6) is a generalized eigenvalue equation. Only if the 
basis of trial wavefunctions is chosen as orthonormal (Eq. 1.10) will S = 1, and the equation 
be an ordinary eigenvalue equation. 

The standard textbook problem of chemical bonding is the formation of molecular orbitals 
from atomic orbitals in describing the molecular ion H}. 

Exercise 12 1-d 

For a potential V(x) = —S(x) — 5{x + a), use <fi mo i(x) = cy4>e{x) + C2<Pe{x + a) as a trial 
wavefunction, and calculate the bonding and antibonding energies as a function of atomic 
separation a. 

3.3 Three dimensions 

There are only slight complications in the real three-dimensional world. We will consider 
only spherical 3-d potentials, v(r), where r = |r|. Then the spatial parts of the orbitals are 
characterized by 3 quantum numbers: 

<f>ni m (r) = Rrd(r)Y kn (6,<p), (3.9) 

where the Y[ m 's are spherical harmonics. Inserting (3.9) into the Schrodinger Equation gives 
the radial equation 

+ ^ +v{r) } Rnl{r) = £ - Al{r) - (3 - 10) 

These are written about in almost all textbooks on quantum mechanics, with emphasis on the 
Coulomb problem, v(r) = —Z/r. Less frequently treated is the three-dimensional harmonic 
oscillator problem, v(r) = kr 2 /2. 

Exercise 13 Three dimensions 

Use the trial wavefunctions of an exponential and a Gaussian in r to deduce the ground-state 
energy and orbital of (a) the hydrogen atom (Z = 1), and (b) the harmonic oscillator with 

k = l. 
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3.4 Differential equations 

We have just seen how to construct various trial wavefunctions and find their parameters by 
minimizing the energy. In a few lucky cases, we've found the exact solution to the given 
problem. 

But in the real world, where analytic solutions are few and far between, we cannot be 
relying on lucky guesses. When we're given a variational problem to solve, we must find the 
exact answer, and be sure we've done so. 

Happily, we can do so, by applying the techniques of Chapter 2 to the variational problem of 
Eq. (3.1). We choose the wavefunction real, and want to minimize E[(j>], given the constraint 
that the wavefunction be normalized. Thus we must minimize: 

F[4>]=T[4>]+V[<t>]-ii{N[$\-l) (3.11) 

where N[(j>] = f dx(j) 2 (x), and p is a Lagrange multiplier. Taking functional derivatives, we 
find: 

-4>"{x), -^— = 2v{x)(j ) {x), -^ = 2<j>(x) (3.12) 



Thus 



5(j)(x) 8(f>(x) Scj)(x) 

SF 



-(j)"{x) + 2(v(x)-p)(f>(x) (3.13) 



5cj>(x) 

At the minimum, this should vanish, yielding: 

- l -f(x)+v{x)0(x)=p<b[x) (3.14) 

and we can identify p = E by integrating both sides with <t>(x). 

This should all make a lot of sense. From the Rayleigh-Ritz variational principle, we have 
rededuced the Schrodinger equation. Note however that all eigenstates satisfy the Schrodinger 
equation, but only the ground state satisfies the variatonal principle (all others are only local 
minima). 

The general rule is that, given a variational principle, trial wavefunctions can be used to 
get approximate answers, but the exact solution can often be found by solving the differential 
equation that results from extremizing the functional, given the constraints. 

3.5 Virial theorem 

We close this chapter with a statement, without proof, of the virial theorem in one dimension. 
The proof will be given later in the book. The virial theorem says that, for an eigenstate of 
a given potential, 

2T=<*£>. (3.15) 
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In the special case in which v(x) = Ax p , (x 



dv 

(h: 



) = p(V), so that 



2T = pV 



(3.16) 



This can save a lot of time when doing calculations, because one only needs either the kinetic 
or the potential energy to deduce the total energy. Note that, for these purposes, the S-fn 
acts as if it had power p = —I. 

The theorem is also true for three-dimension problems, and also true for interacting prob- 
lems, once the potential is interpreted as all the potential energy in the problem. Thus, for 
any atom, p = — 1, and E = —T = —V/2, where V is the sum of the external potential and 
the electron-electron repulsion. 

Exercise 14 Virial theorem for one electron 

Check Eq. (3.16) for the ground state of (a) the Id H atom, (b) the Id harmonic oscillator, 
(c) the three-dimensional Hydrogen atom. Also check it for Ex. (9) and (10). 

Some approximations automatically satisfy the virial theorem, others do not. In cases that 
do not, we can use the virial theorem to judge how accurate the solution is. In cases that do, 
we can use the virial theorem to check we dont have an error, or to avoid some work. 

3.6 Questions 

All the questions below are conceptual. 

1. Suggest a good trial wavefunction for a potential that consists of a negative delta function 
in the middle of a box of width L. 

2. Sketch, on the same figure, the Id H-atom for Z = 1/2,1, and 2. What happens as 

Z — > oc and as Z — > 0? 

3. What is the exact kinetic energy density functional for one electron in one-dimension? 

4. How would you check to see if your trial wavefunction is the exact ground-state wave- 
function for your problem? 

5. Consider what happens to the potential v(x) = Ax p as p gets large. Can the virial 
theorem be applied to the particle in the box (even if it yields trivial answers)? 
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Two electrons 



1 In this chapter, we review the traditional wavefunction picture of Schrddinger for two 
electrons. We keep everything as elementary as possible, avoiding sophistry such as the 
interaction picture, second quantization, Matsubara Green's functions, etc. These are all 
valuable tools for studying advanced quantum mechanics, but are unneccessary for the 
basic logic of our presentation. 

4.1 Antisymmetry and spin 

When there is more than one particle in the system, the Hamiltonian and the wavefunction 
includes one coordinate for each particle, as in Eqn. (1.1) for the H 2 molecule. Furthermore, 
electrons have two possible spin states, up or down, and so the wavefunction is a function both 
of spatial coordinates and spin coordinates. A general principle of many-electron quantum 
mechanics is that the wavefunction must be antisymmetric under interchange of any two sets 
of coordinates. All this is covered in elementary textbooks. For our purposes, this implies the 
2-electron wavefunction satisfies: 

#(r 2 , cr 2 , ri, <7i) = -*(r b <r u r 2 , er 2 ). (4.1) 

The ground-states of our systems will be spin-singlets, meaning the spin-part of their wave- 
function will be antisymmetric, while their spatial part is symmetric: 

(Fin, r 2 a 2 ) = $(n, r 2 ) X Singlet (^, <r 2 ), (4.2) 

where the spatial part is symmetric under exchange of spatial coordinates and the spin part 

X Singlet (o- 1 ,a 2 ) = ^= 2 (\ Tl )- | IT )) is antisymmetric. 

4.2 Hartree-Fock 

For more than one electron, the operators in the Hamiltonian depend on all the coordinates, 
but in a simple way. Both the kinetic energy and the external potential are one-body operators, 

1 ©2000 by Kieron Burke. All rights reserved. 
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meaning that they are sums of terms which each depend on only one coordinate at a time: 

1 ( d 2 d 2 \ 

T = ~2 I dxf + dxjj Vext = Vext ^ + v ^ x ^ ( 4 - 3 ) 

The electron-electron repulsion operator is a two-body operator, each term depending on 
two coordinates simultaneously. Note that it is this term that complicates the problem, 
by coupling the two coordinates together. The interaction between two electrons in three 
dimensions is Coulombic, i.e., l/|r — r'|. This is homogeneous of degree -1 in coordinate 
scaling. However, l/\x — x'\ in one-dimension is an exceedingly strong attraction, and we 
prefer to use 5(x — x'), which has the same scaling property, but is much weaker and more 
short-ranged. 

We will now consider the problem of one-dimensional He as a prototype for two electron 
problems in general. The Hamiltionian is then: 

H = -\ (J^ + |||) - Z6( X1 ) - Z5(x 2 ) + S( X1 - x 2 ). (4.4) 

To find an exact solution to this, we might want to solve the Schrodinger equation with 
this Hamiltonian, but to find approximate solutions, its much easier to use the variational 
principle: 

E = mm(q\H\t>), J dx Y dx 2 \^{x u x 2 ) \ 2 = 1. (4.5) 

The most naive approach might be to ignore the electron-electron repulsion altogether. 
Then the differential equation decouples, and we can write: 

V(x 1 ,x 2 ) = 4>{xi)<j>{x 2 ) (ignoring l/ ee ) (4.6) 

where both orbitals are the same, and satisfy the one-electron problem. We know from 
chapter 3 that this yields 

<t>(x) = y/Zexp{-Z\x\). (4.7) 

This decoupling of the coordinates makes it possible to handle very large systems, since we 
need only solve for one electron at a time. However, we have made a very crude approximation 
to do this. The contribution of the kinetic energy and the potential energy is then just twice 
as big as in the single electron system: 

T, = 2(Z 2 /2) V ext = 2(-Z 2 ) V ee = (4.8) 

and therefore the total energy becomes 

E = 2e = -Z 2 (4.9) 

which gives for He (Z=2) E = —4. This is in fact fewer than the true ground-state energy 
of this problem, because we have failed to evaluate part of the energy in the Hamiltonian. 
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To improve on our estimate, and restore the variational principle, we should evaluate the 
expectation value of V ee on our wavefunction: 

V ee = r dx r dx'Wx)\ 2 \4>{x')\ 2 5(x - x') = Z 2 f°° dx exp(-4Z|x|) = % . (4.10) 

J— oo J— oo J— oo 2 

For Helium (Z=2), V ee = +1 so that the ground state energy becomes E = — 4 + 1 = —3. 
Because we evaluated all parts of the Hamiltonian on our trial wavefunction, we know that 
the true E < -3. 

However, using the variational principle, we see that, with almost no extra work, we can 
do better still. We chose the length scale of our orbital to be that of the V ee = problem, 
i.e., a = Z, the nuclear charge. But we can treat this instead as an adjustable parameter, 
and ask what value minimizes the energy. Inserting this orbital into the components of the 
energy, we find: 

T, = 2(a 2 /2) V ext = -2Za V ee = a/2 (4.11) 
so that the total energy, as a function of a, is 

E(a) = a 2 - 2a(Z- 1/4). (4.12) 
Minimizing this, we find a min = Z — 1/4 and thus 

Emin = -a 2 min = -(Z - 1/4) 2 = - Q 2 = -3.0625. (4.13) 

We have lowered the energy by 1/16 of a Hartree, which may not seem like much, but is 
about 1.7 eV or 40 kcal/mol. Chemical accuracy requires errors of about 1 or 2 kcal/mol. 

The best solution to this problem is to find the orbital that produces the lowest energy. 
We could do this by including many variational parameters in a trial orbital, and minimize 
the energy with respect to each of them, giving up when addition of further parameters 
has negligible effect. A systematic approach to this would be to consider an infinite set 
of functions, usually of increasing kinetic energy, and to include more and more of them. 
However, having learned some functional calculus, we can provide a more direct scheme. 

The energy, as a functional of the orbital, is: 

E[4>\ = 2T S [0] + 2V ext [cf\ + U[4>]/2, U[4>\ = \r dxn 2 (x) (4.14) 

Z ■> — oo 

where n(x) = 2\<p(x)\ 2 , T s and V ex t are the one-electron functionals mentioned in chapter 3, 
and U is the 1-d equivalent of the Hartree energy (discussed more below). For two electrons 
in HF, V ee = U/2, and we discuss this much more later. If we simply minimize this functional 
of the orbital, subject to the restriction that the orbital is normalized, we find: 

~H ~ zs{x} + l vu{x) } = e0(x) ' (4 - i5) 
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where 

vJx) = = r dx' n(x')5(x - x') = n(x) (4.16) 

on[x) " , -° c 

is the Hartree potential. This is also known as the classical or electrostatic potential, as it 
is the electrostatic potential due to the charge distribution in classical electrostatics (if the 
electrostatic interaction were a 5-function, for this case). (Compare with Eq. (1.6) of the 
introduction.) This equation is the Hartree-Fock equation for this problem. 

Exercise 15 Hartree-Fock equations for two electrons 

Derive the Hartree-Fock equation for two electrons, Eq. (4.15), by minimizing the energy as 
a functional of the orbital, Eq. (4.14), keeping the orbital normalized, using the techniques 
of Chapter ??. 

There are several important aspects of this equation that require comment. First, note the 
factor of 2 dividing the Hartree potential in the equation. This is due to exchange effects: 
The Hartree potential is the electrostatic potential generated by the total charge density due 
to all the electrons. A single electron should only see the potential due to the other electron, 
hence the division by two. 

Next, note that these are self-consistent equations, since the potential depends on the 
density, which in turn depends on the orbital, which is the solution of the equation, etc. 
These can be solved in practice by the iterative method described in the introduction. 

To get an idea of what the effective potential looks like, we construct it for our approximate 
solution above. There a = 1.75 for He, so the potential is a delta function with an added 
Hartree potential of 4exp(— 7\x\). This contribution is positive, pushing the electron farther 
away from the nucleus. We say the other electron is screening the nucleus. The effective 
nuclear charge is no longer Z, but rather Z — 1/4. 

It is straightforward to numerically solve the orbital equation Eq. (4.15) on a grid, essen- 
tially exactly. But there is an underlying analytic solution in this case. Write 

= 7 cosech (7 \x\ + /?) (4.17) 

where coth(S) = Z/j, and 7 = Z — 1/2. Insertion into Eq. (4.15) yields 

e = - 7 2 /2, (4.18) 

but this does not mean E = — 7 2 . The energy of the effective non-interacting system is not 
the energy of the original interacting system. In this case, from the derivation, we see E = 
2e-U/2 (since e = T s + V mt + U/2). We find U/2 = Z/2-1/Q, or E = -Z 2 + Z/2-l/\2 as 
the exact HF energy, slightly lower than our crude estimate using hydrogenic wavefunctions. 
But this slight difference is 1/48, or 0.6 eV or 13 kcal/mol, which is still a significant error 
by modern standards. 

In Fig. 4.1, we plot both the exact orbital and the scaled Hydrogenic orbital, showing there 
is very little difference. 




0.5 1 1.5 2 2.5 3 3.5 4 

X 

Figure 4.1: HF orbital for Id He atom, both exact and simple exponential with effective charge. 



Exercise 16 Analytic solution to HF for Id He 

This exercise is best done with a table of integrals, a mathematical symbolic program, or 
simply numerically on a grid. 

1. Show that the orbital of Eq. (4.17) is normalized. 

2. Show that the orbital of Eq. (4.17) satisfies the HF equation with the correct eigenvalue. 

3. Check that the HF solution satisfies the virial theorem ( for power -1 ). 

Exercise 17 Ionization energy of 1-d He 
The ionization energy of He is defined as 

I = E X -E 2 (4.19) 

where E N is the ground-state energy with N electrons in the potential. 

1. What is the estimate of the ionization energy given by the crude Hydrogenic orbital 
approximation? 

2. What is the ionization energy given by the HF solution? 

3. Can you think of another way to estimate the ionization energy from the HF solution? 
Exercise 18 He in 3-d 

Repeat the approximate HF calculation for He in three dimensions. What is the orbital energy 
the effective nuclear charge, and the total physical energy? Make a rigorous statement about 
both the exact HF energy and the true ground-state energy. Plot the effective potential 
(external plus half Hartree) using your approximate solution. (For help, this calculation is 
done in many basic quantum books). 
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atom 




E 




E c 


H- 


-0.486 


-0.528 


-0.381 


-0.042 


He 


-2.862 


-2.904 


-1.025 


-0.042 


Be++ 


-13.612 


-13.656 


-2.277 


-0.044 



Table 4.1: Energies for two-electron ions. 

4.3 Correlation 

The approximate solution of the Hartree-Fock equations is not exact, because the true wave- 
function is not a product of two orbitals, but is rather a complicated function of both vari- 
ables simultaneously. The true wavefunction satisfies the exact Schrodinger equation, and 
also minimizes the ground-state energy functional for the given external potential. In tra- 
ditional quantum chemistry, the correlation energy is defined as the difference between the 
Hartree-Fock energy and the exact ground-state energy, i.e., 

£trad = E _ E m ( 4 20 ) 

We illustrate the effects of correlation on the simplest system we have, He. We can 
adjust the relative importance of correlation by considering various values of Z . Perhaps 
counterintuitively, for large Z, the external potential dominates over the repulsion, and so the 
non-interacting solution becomes highly accurate. As Z is reduced, the interaction becomes 
more important. For larger Z values, Hartree and exchange dominate over correlation. But 
when Z is of order 1, correlation becomes important. Eventually, if Z is made too small, the 
system becomes unstable to ejecting one of the electrons, i.e., its energy becomes less than 
that of the one-electron ion, i.e., its ionization potential passes through zero. For Z just a 
little greater than that, the system has significant correlation. We illustrate these results in 
Table 4.3. For large Z, the total energy grows as — Z 2 (the Hydrogenic result), while the 
exchange energy grows as — 5Z/8. Thus the magnitude of the exchange energy grows with 
Z, but its relative size diminishes. Similarly, clearly E c is almost independent of Z. So E c 
becomes an ever smaller fraction of E x as Z grows. Finally, notice that, for H~ , the exchange 
energy is so small that correlation is 11% of it. The two-electron ion unbinds at the critical 
value of Z = 0.911. 

Returning to the Id world, all the same remarks apply qualitatively. The ground-state 
energy for Id He is -3.154. Since the HF energy is exactly -3 1/12, the correlation energy is 
71 mH. Thus we need to estimate correlation energies to within about 10% for useful accuracy 
in quantum chemistry. This may seem daunting, but we've made great strides so far: Even 
our crudest method was good to within 5%. This is a gift of the variational principle, which 
implies that all errors in the ground-state energy are at least second-order in the error in the 
wavefunction (don't expect other expectation values to be as good). Its just that we need 
about 0.2%! 
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fixed Z adjustable Z 



HF exact 



E 



-3 -3.0625 



-3.0833 -3.154 



% err 



-5 -3 



-2 



Table 4.2: Approximate Energies for Id He. 



Exercise 19 Ions on the verge of ionization 

1. Using the simple exponential approximation, calculate the critical Z at which Id He 
becomes unstable. 

2. Repeat using the exact HF result. Is there an alternative method of estimating Z c in this 
case? Which is better? 

3. The exact energy is accurately reproduced by 



where a = 3/4 - 4/(3tt) = 0.3255868 and (5 = 0.0342. Find the critical value of Z. 
4. Calculate the ratio of E c to E x at Z c . 

4.4 Questions 

All the questions below are conceptual. 

1. Is the HF estimate of the ionization potential for 1-d He an overestimate or underesti- 
mate? 

2. Consider the approximate HF calculation given in section 4.2. Comment on what it does 
right, and what it does wrong. Suggest a simple improvement. 

3. Which is bigger, the kinetic energy of the true wavefunction or that of the HF wavefunc- 
tion for Id He? 




(4.21) 



Chapter 5 
Many electrons 



In this chapter, we generalize the concepts and ideas introduced for two electrons to the 
N electron case. 

5.1 Ground state 

We define carefully our notation for electronic systems with N electrons. First, we note 
that the wavefunction for TV electrons is a function of 37V spatial coordinates and TV spin 
coordinates. Writing x t = (r,,cr,;) to incorporate both, we normalize our wavefunction by: 

j dx\ . ..J dx N \V(xi, . . . ,x N )\ 2 = 1, (5.1) 

where / dx denotes the integral over all space and sum over both spins. Note that the 
antisymmetry principle implies 

. . .,Xj, ...,x h ...) = .,Xj, . . .). (5.2) 

The electronic density is defined by 

n(r) =iVE fdx 2 ... f dx N |*(r, a, x 2 , ■ ■ ■ , x N )\ 2 , (5.3) 

and retains the interpretation that n(r)d 3 r is the probability density for finding any electron 
in a region d 3 r around r. The density is normalized to the number of electrons 

j d 3 r n(r) = N. (5.4) 

Our favorite operators become sums over one- and two-particle operators: 

1 N 

T = -^l (5-5) 

1 1=1 

V ext = E u ext(ri), (5.6) 

i=\ 

and 

1 N 1 

Ke = 2 E rr— n (5 - 7) 

47 
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Note that the factor of 2 in the electron-electron repulsion is due to the sum running over all 
pairs, e.g., including (1,2) and (2,1). The Schrodinger equation is then 

(f + t/ ee + t/ ext ) *{x 1 ,...,xit) = E9(x 1 ,...,x N ). (5.8) 

The ground-state energy can be extracted from the variational principle: 

75 = mm(*|f + V; e + V; xt |*), (5.9) 

once the minimization is performed over all normalized antisymmetic wavefunctions. 



5.2 Hartree-Fock 

In the special case of non-interacting particles, we denote the wavefunction by $ instead of 
^, and this will usually be a single Slater determinant of occupied orbitals, i.e., 



$(xi, . . . ,xn) 



4>\{x N ) ■■■ 4> N (x N ) 



(5.10) 



For systems with equal numbers of up and down particles in a spin-independent external 
potential, the full orbitals can be written as a product, (j>i(x) = <j>i{r)\a), and each spatial 
orbital appears twice. The Hartree-Fock energy is 

E = mm(^\f + V ee + V ext \^}, (5.11) 

The first contribution is just the non-interacting kinetic energy of the many orbitals, and the 
last is their external potential energy. 

The Coulomb interaction for a single Slater determinant yields, due to the antisymmetric 
nature of the determinant, two contributions: 

<$|\> ee |$> = [/[$] + £ x [$]. (5.12) 

The first of these is called the direct or Coulomb or electrostatic or classical or Hartree 
contribution, given in Eq. F.4, with the density being the sum of the squares of the occupied 
orbitals. This is the electrostatic energy of the charge density in electromagnetic theory, 
ignoring its quantum origin. The second is the Fock or exchange integral, being 

em =-ki f d *r j# r . m*m yw) (5 . 13) 

for a determinant of doubly-occupied orbitals. This is a purely Pauli-exclusion principle effect. 
Exercise 20 Exchange energies for one and two electrons 

Argue that, for one electron, E x = —U, and show that, for two electrons in a singlet, 

£x = -U/2. 
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-0.5 
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-37.689 


-37.845 


-0.156 
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8 


-74.809 


-75.067 


-0.258 


F 


9 


-99.409 


-99.733 


-0.324 


Ne 


10 


-128.547 


-128.937 


-0.39 


Ar 


18 


-526.817 


-527.539 


-0.722 


Kr 


36 


-2752.055 


-2753.94 


-1.89 


Xe 


54 


-7232.138 


-7235.23 


-3.09 


Rn 


86 


-22866.745 


-22872.5 


-5.74 



Table 5.1: Total energies in HF, exactly, and correlation energies across first row and for noble gas atoms. 

To find the orbitals in Hartree-Fock, we must minimize the energy as a functional of 
each orbital <j>i{r), subject to the constraint of orthonormal orbitals. Doing this, using the 
techniques of chapter 2 yields the Hartree-Fock equations 

(-^V 2 + v ext (r) + w H (r)) &(r) + /«[{4}](r) = £j <A(r) (5.14) 
where the last contribution to the potential is due to the Fock operator, and is defined by 

f F .mm = -£/ dv ^y fr(r). (5.15) 

Note that this odd-looking animal is or6ito/-dependent, i.e., it is a different function of r for 
each occupied orbital. Also, for just one electron, the Hartree and Fock terms in the potential 
cancel, as they should. In order to discuss the pro's and con's of HF calculations, we must 
first discuss the errors it makes. 

5.3 Correlation 

The definition of correlation energy remains the same for N electrons as for two: It is the 
error made by a Hartree-Fock calculation. Table 5.3 lists a few correlation energies for 
atoms. We see that correlation energies are a very small (but utterly vital) fraction of the 
total energy of systems. They are usually about 20-40 mH/electron, a result we will derive 
later. 
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Exercise 21 Atomic energies 

This exercise uses the numbers in Table 5.3. 

1. By plotting \n(—E) versus ln(Z), find the dependence of the total energy on Z for large 
Z. 

2. Repeat for the correlation energy. 

3. Plot the correlation energy (a) across the first row and (b) down the last column of the 
periodic table. Comment on the results. 

But we strive for chemical accuracy in our approximate solutions, i.e., errors of less than 
2 mH per bond. Another important point to note is that in fact, we often do not need total 
energies to this level of accuracy, but rather only energy differences, e.g., between a molecule 
and its separated constituent atoms, in which the correlation contribution might be a much 
larger fraction, and in which compensating errors might occur in the separate calculation of 
the molecule and the atoms. An example of this is the core electrons, those electrons in closed 
subshells with energies below the valence electrons, which are often relatively unchanged in 
a chemical reaction. A large energy error in their contribution is irrelevant, as it cancels out 
of energy diffferences. 

Quantum chemistry has developed many interesting ways in which to calculate the corre- 
lation energy. These can mostly be divided into two major types: perturbation theory, and 
wavefunction calculations. The first is usually in the form of Moller-Plesset perturbation the- 
ory, which treats the Hartree-Fock solution as the starting point, and performs perturbation 
theory in the Coulomb interaction. These calculations are non-variational, and may produce 
energies below the ground-state energy. Perhaps the most common type of wavefunction 
calculation is configuration interaction (CI). A trial wavefunction is formed as a linear com- 
bintation of products of HF orbitals, including excited orbitals, and the energy minimized. 
In electronic structure calculations of weakly-correlated solids, most often density functional 
methods are used. Green's function methods are often applied to strongly-correlated systems. 

Exercise 22 Correlation energy 

Show that the correlation energy is never positive. When is it zero? 

An interesting paradox to note in chemistry is that most modern chemists think of reac- 
tivity in terms of frontier orbitals, i.e., the HOMO (highest occupied molecular orbital) and 
the LUMO (lowest unoccupied), and their energetic separation, as in Fig. 1.3. In the wave- 
function approach, these are entirely constructs of the HF approximation, which makes errors 
of typically 0.2 Hartrees in binding energies. Thus, although this approximation obviously 
contains basic chemical information, needed for insight into chemical reactivity, the resulting 
thermochemistry is pretty bad. On the other hand, to obtain better energetics, one adds 
many more terms to the approximate wavefunction, losing the orbital description completely. 
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We will see later how DFT resolves this paradox, by showing how an orbital calculation can 
in principle yield the exact energetics. 

5.4 Atomic configurations 

We have seen how the Hartree-Fock equations produce a self-consistent set of orbitals with 
orbital energies. In many cases, the true wavefunction has a strong overlap with the HF 
wavefunction, and so the behavior of the interacting system can be understood in terms of 
the HF orbitals. We occupy the orbitals in order. 

This is especially true for the atoms, and the positions of the HF orbitals largely determine 
the strucutre of the periodic table. If we first ignore interaction, and consider the hydrogenic 
levels, we get the basic idea. For Hydrogenic atoms (N = 1), there is an exact degeneracy 
between all orbitals of a given principal quantum number. For each n, there are n values of 
I, ranging from to n — 1. For each I, there are 21 + 1 values of m; = — I, 0, ..I. Thus 
each closed shell contains 

n-1 

9n = £(2/ + 1) = n(n- 1) +n = n 2 (5.16) 

1=0 

Since each orbital can hold 2 electrons, the nth-shell contains 2n 2 electrons. The first shell 
closes at N = 2, the next at 10, the next at 28, and so forth. 

This scheme works up to Ar, but then fails to account for the transition metals. To 
understand these, we must consider the actual HF orbitals, which are not hydrogenic. Since 
the HF potentials are not Coulombic, the degeneracy w.r.t. different /-values is lifted, and 
we have e„i. As expected, for fixed n, the higher I, the less negative the energy. Thus the 2s 
orbitals lie lower than 2p. This effect does not change the filling order described above. 

But for Z between 7 and 20, the 4s orbital dips below the 3d. So it gets filled before the 
3d, and, for example, the ground-state configuration of Ca is [Ar]4s 2 . Then the 3d starts 
filling, and so the transition metals appear a little late. Even when the 3d orbital energy drops 
below the 4s, the total energy remains lower with the 4s filled. The situation is similar for 
higher n. 

A list of the filled orbitals is called the electronic configuration of a system. This still does 
not determine the ground-state entirely, as many different angular momenta combinations, 
both spin and orbital, (terms) can have the same configuration. Hund's rules are used to 
choose which one has the lowest energy. Our treatment will not require details beyond 
knowing the lowest configuration. 

5.5 Atomic densities 

Since we will be constructing a formally exact theory of interacting quantum mechanics based 
on the one-electron density, it seems appropriate here to discuss and show some pictures of 
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this quantity for atoms. All atomic densities, when plotted as function of r with the nucleus 

25 




Figure 5.1: Radial density in Ar atom. 

at the origin, have been found to decay monotonically, although that's never been generally 
proven. Close to a nucleus, the external potential dominates in the Schrbdinger equation, 
and causes a cusp in the density, whose size is proportional to the nuclear charge: 



-2Z n(0) 



(5.17) 



dr lr = ' 

for a nucleus of charge Z at the origin. This is Kato's cusp condition, and we already saw 
examples of this in earlier exercises. It is also true in our one-dimensional examples with 
delta-function external potentials. 

Typically, spherical densities are plotted multiplied by the phase-space factor Airr 2 . This 
means the area under the curve is precisely N. Furthermore, the different electronic shells 
are easily visible. We take the Ar atom as an example. In Fig. 5.1, we plot its radial density. 
This integrates to 18 electrons. We find that the integral up to r = 0.13 contains 2 electrons, 
up to 0.25 contains 4, 0.722 contains 10, and 1.13 contains 12. These correspond to the 2 
Is electrons, 2 2s electrons, 6 2p electrons, and 2 3s electrons, respectively. Note that the 
peaks and dips in the radial density roughly correspond to these shells. 

The decay at large distances is far more interesting. When one coordinate in a wavefunction 
is taken to large distances from the nuclei, the iV-electron ground-state wavefunction collapses 
to the product of the square-root of the density times the (N— l)-electron wavefunction. This 
means the square-root of the density satisfies a Schrodinger-like equation, whose eigenvalue 
is the difference in energies between the two systems: 

\Jn(r) = Ar f) exp(— ar) (r — > oo) (5.18) 

where a = \f2I, and 

/ = E(N — 1) — E(N) (5.19) 
is first ionization potential. In fact, the power can also be deduced, (3 = (Z — N + I) /a — 1, 
and A is some constant. Thus, a useful and sensitive function of the density to plot for 
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spherical systems is 

«M = (i/2)^W), (520) 

since k(0) = —Z, while /-c(r) — > —a, as r — > oo. The ground-state density can always be 
reconstructed from 

n(r) = Cexp(jf°° dr'2 K (r')) (5.21) 

where the constant is determined by normalization. 

In a one-dimensional world with delta function interactions, these conditions remain true, 
except for the details of the power law in front of the exponential in Eq. (5.18). In particular, 
we have already seen the cusp in the orbital for 1-d hydrogenic atoms. We may see dramatic 
evidence for these conditions in our HF solution of the 1-d He atom, by plotting k(x) on the 
left of Fig. 5.2. Near the nucleus, the exact cusp condition is satisfied, so that k equals -2, 
but at large distances, it tends to a smaller constant, determined by the ionization potential. 
We can also clearly see the difference between the approximate and exact HF solutions here, 
which was not so visible in the last figure. On the right, we have also plotted n(r) for real 
He, and see that the effect of correlation on the density is extremely small: The HF and exact 
densities are very close. 
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Figure 5.2: k in (a) Id He atom in HF approximation, both exact and approximate; (b) in real He, in HF and exactly 

5.6 Questions 

All the questions below are conceptual. 

1. Does each orbital in a HF calculation satisfy the same Schrodinger equation? 

2. What is the relationship between E UF and E-Ii ej 

3. Speculate on why the correlation energy of Li is about the same as that of He. 
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4. Which is bigger, the kinetic energy of the true wavefunction or that of the HF wavefunc- 
tion, for an atom? 

5. What do you expect happens to the correlation energy of the two-electron ions as Z — > 

oc? 

6. What do you expect happens to the correlation energy of the four-electron ions as Z — > 

oc? 



Part II 
Basics 
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Chapter 6 

Density functional theory 



1 This chapter deals with the foundation of modern density functional theory as an exact 
approach (in principle) to systems of interacting particles. In our case, electrons, i.e., 
fermions interacting via the Coulomb repulsion. 

6.1 One electron 

To illustrate the basic idea behind density functional theory, we first formulate the problem 
of one electron as a density functional problem. We know that the ground-state energy can 
be written as a density functional instead of as a wavefunction functional. This functional is 

E[n] = 1 f d 3 r + / d 3 rv(r) n(r) (6.1) 

in any number of dimensions. We must minimize it, subject to the constraint that / d 3 r n(r) = 
1. Using Lagrange multipliers, we construct the auxiliary functional: 

H[n] = E[n]-fiN (6.2) 
where N = J d i rn(r). Minimizing this yields 

SH ST VW , N „ 

oW) = Hr) +V[r) -" = ° (6 - 3) 
Using the derivative from Ex. 5, we find the Schrodinger equation for the density: 

+ + «(r))n(r)=/m(r) (6.4) 

with boundary conditions that n(r) and Vn(r) — > at the edges. We can identify the 
Lagrange multiplier by integrating both sides over all space at the solution. Since the integral 
of any Laplacian vanishes (with the given boundary conditions), we see that p, = E. 

Exercise 23 Checking equation for density 

Show that the ground-state density for a particle in a box satisfies Eq. (6.4). 

1 ©2000 by Kieron Burke. All rights reserved. 
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6.2 Hohenberg-Kohn theorems 

It is self-evident that the external potential in principle determines all the properties of the sys- 
tem: this is the normal approach to quantum mechanical problems, by solving the Schrodinger 
equation for the eigenstates of the system. 

The first Hohenberg-Kohn theorem demonstrates that the density may be used in place 
of the potential as the basic function uniquely characterizing the system. It may be stated 
as: the ground-state density n(r) uniquely determines the potential, up to an arbitrary 
constant. 

In the original Hohenberg-Kohn paper, this theorem is proven for densities with non- 
degenerate ground states. The proof is elementary, and by contradiction. Suppose there 
existed two potentials differing by more than a constant, yielding the same density. These 
would have two different ground-state wavefunctions, and \t 2 - Consider 3' 2 as a trial 
wavefunction for potential v ext .i(r). Then, by the variational principle, 

(* 2 |T+V; e + V; xt;1 |* 2 ) > (*l|T+Vee + Vert,i|*i>. (6.5) 

But since both wavefunctions have the same density, this implies 

<* 2 |T + f ee |* 2 ) > {^If+Veel*!}- (6-6) 

But we can always swap which wavefunction we call 1 and which we call 2, which reverses 
this inequality, leading to a contradiction, unless the total energies of the two wavefunctions 
are the same, which implies they are the same wavefunction by the variational principle and 
the assumption of non-degeneracy. Then, simple inversion of the Schrodinger equation, just 
as we have done several times for the one-electron case, yields 

N 1 N 11 

Ea^-.EV^+.E^- 1 (6.7) 

i=l Z i=l L j^i \ v i — r j 

This determines the potential up to a constant. 

Exercise 24 Finding the potential from the density 

For the smoothed exponential, <p(x) = C(l + \x\) exp(— \x\), find the potential for which this 
is an eigenstate, and plot it. Is it the ground state of this potential? 

An elegant constructive proof was found later by Levy, which automatically includes degen- 
erate states. It is an example of the constrained search formalism. Consider all wavefunctions 
which yield a certain density r?,(r). Define the functional 

F[n] = min(* If + V ee \ *) (6.8) 

where the search is over all antisymmetric wavefunctions yielding r?,(r). Then, for any n(r), 
any wavefunction minimizing T + V ee is a ground-state wavefunction, since the ground-state 
energy is simply 

E = mm (F[n] + / d 3 r v ext (r) n(r)) , (6.9) 
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from the variational principle, where the search is over all normalized positive densities. Any 
such wavefunction can then be fed into Eq. (6.7) to construct the unique corresponding 
potential. We denote the minimizing wavefunction in Eq. (6.8) by ty[n]. This gives us a 
verbal definition of the ground-state wavefunction. The exact ground-state wavefunction of 
density n(r) is that wavefunction that yields n(r) and has minimizes T + V ee . We may 
also define the exact kinetic energy functional as 

T[n] = (tf[ra]|fjtt[n]), (6.10) 

and the exact electron-electron repulsion functional as 

V ee [n] = (y[n]\v ee \v[n]}. (6.11) 

The second Hohenberg-Kohn theorem states that the functional F[n] is universal, i.e., 
it is the same functional for all electronic structure problems. This is evident from Eq. 
(6.8), which contains no mention of the external potential. To understand the content of 
this statement, consider some of our previous problems. Recall our orbital treatment of 
one-electron problems. The kinetic energy functional, T[</>] = \\ dx \<p'(x)\ 2 , is the same 
functional for all one-electron problems. When we evaluate the kinetic energy for a given 
trial orbital, it is the same for that orbital, regardless of the particular problem being solved. 
Similarly, when we treated the two electron case within the Hartree-Fock approximation, we 
approximated 

F[n] * F H >] = 1 / d'VM + I / rfV / rfV H< : r)n{r '\ (6.12) 
8 J n 4 J J \r — r '| 

Then minimization of the energy functional in Eq. (6.9) with F HF [n] yields precisely the 
Hartree-Fock equations for the two-electron problem. The important point is that this single 
formula for F HF [n] is all that is needed for any two-electron Hartree-Fock problem. The 
second Hohenberg-Kohn theorem states that there is a single F[n] which is exact for all 
electronic problems. 

Exercise 25 Errors in Hartree-Fock functional 

Comment, as fully as you can, on the errors in F HF [n] relative to the exact F[n] for two 
electrons. 

The last part of the Hohenberg-Kohn theorem is the Euler-Lagrange equation for the 
energy. We wish to mininize E[n] for a given v e xt(r) keeping the particle number fixed. We 
therefore minimize E[n] — fiN, as in chapter ??, and find the Euler-Lagrange equation: 

8F 



We can identify the constant /i as the chemical potential of the system, since ji = dE/dN. 
The exact density is such that it makes the functional derivative of F exactly equal to the 
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negative of the external potential (up to a constant). Note that it would be marvellous if 
we could find an adequate approximation to F for our purposes, so that we could solve Eq. 
(6.13) directly. It would yield a single integrodifferential equation to be solved, probably 
by a self-consistent procedure, for the density, which could then be normalized and inserted 
back into the functional E[n], to recover the ground-state energy. In the next section, we 
will examine the original crude attempt to do this (Thomas-Fermi theory), and find that, 
although the overall trends are sound, the accuracy is insufficient for modern chemistry and 
materials science. Note also that insertion of F hf [ti] will yield an equation for the density 
equivalent to the orbital HF equation. 

An important formal question is that of t;-representability. The original HK theorem was 
proven only for densities that were ground-state densities of some interacting electronic prob- 
lem. The constrained search formulation extended this to any densities extracted from a 
single wavefunction, and this domain has been further extended to densities which result 
from ensembles of wavefunctions. The interested reader is referred to Dreizler and Gross for 
a thorough discussion. 

6.3 Thomas-Fermi 

In this section, we briefly discuss the first density functional theory (1927), its successes, and 
its limitations. Note that it predates Hartree-Fock by three years. 

In the Thomas-Fermi theory, F[n] is approximated by the local approximation for the 
(non-interacting) kinetic energy of a uniform gas, plus the Hartree energy 

F T >] = As / d^\r) + \jd*r J dV . (6.14) 

Several points need to be clarified. First, these expressions were developed for a spin- 
unpolarized system, i.e., one with equal numbers of up and down spin electrons, in a spin- 
independent external potential. Second, in the kinetic energy, we saw in Chapter ?? how the 
power of n can be deduced by dimensional analysis, while the coefficient is chosen to agree 
with that of a uniform gas, yielding As = (3/10)(37r 2 ) 2,/3 . We will derive these numbers in 
detail in Chapter 8. They can also be derived by classical arguments applied to the electronic 
fluid. 

Insertion of this approximate F into the Euler-Lagrange equation yields the Thomas-Fermi 
equation: 

^^(r) + / dV^- + , ext (r) = ,i (6.15) 
3 | r — r' | 

We will focus on its solution for spherical atoms, which can easily be achieved numerically with 

a simple ordinary differential equation solver. We can see immediately some problems. As 

r — > 0, because of the singular Coulomb external potential, the density is singular, n — ► 1/r 3 / 2 

(although still integrable, because of the phase-space factor 4nr 2 ). Again, as r — > oo, 
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the density decays with a power law, not exponentially. But overall, the trends will be 
approximately right. 

The Thomas-Fermi equation for neutral densities, Eq. (6.15), is usually given in term of a 
dimensionless function satisfying: 



\ 



$3 

X 



where $(0) = 1 and $(oc) = 0. The distance x = Z l l' i, rja, where 

1 (^\ 2 '\ 



and the density is given by 



2 V 4 J 

n(r) 



0.885341. 



z 2 /$\ 3 / 2 



(6.16) 



(6.17) 



4na 3 \x } ^ ^ 

Its been found that at the solution for neutral atoms, $'(0) = —B, where B = 1.58807102261, 
and the asymptotic behavior is the trivial solution: 

144 



Mx) 



.3 ' 



(6.19) 



These results, and the density of Fig. 6.1, are needed for the next exercise. 

90 r 
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r 

Figure 6.1: Radial density in the Xe atom, calculated using LDA. 

Exercise 26 Thomas-Fermi for atoms 



1. Using the differential equation and integrating by parts, find the TF kinetic energy. Using 
the virial theorem for atoms, show that 

E = - 3 - Z~'I' A = -.7845 Z 1 ' 3 (6.20) 
7 a y ' 

This is the exact limiting behavior for non-relativistic atoms as Z — > oc. 
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2. Show that the external potential goes as 

V ext = -- Z 7/i = 1.79374 Z 7ri (6.21) 

3. By using the virial theorem, deduce the behavior of the Hartree energy as Z — > oc. 

4. Using table 5.3, compare with energies along first row and for the noble gas atoms. 
Comment. 

5. An accurate radial density for the Xe atom is given in Fig. 6.1. An accurate parametriza- 
tion of $ is 

<S>(y) = (l + ay + 0y 2 exp(- 7 y)f exp(-2a y) (6.22) 

where y = y/x and a = 0.7280642371,/? = -0.5430794693, and 7 = 0.3612163121. 
Plot the approximate TF density on the same figure, and comment. 

6. Calculate the expectation value (1/r) for Thomas-Fermi. Do atomic radii grow or shrink 
with Z? Explain why. 

To see the quality of our TF solution, we calculate the n{r) function of the last chapter, 
and compare with that of the exact solution. We see that, although it messes up all the 
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Figure 6.2: «(r) for the Ar atom, both exactly and in Thomas-Fermi calculation of Exercise ??. 

details (nuclear cusp, shell structure, decay at large distances), it has a remarkably good 
overall shape for such a simple theory. 

Lastly, Teller showed that molecules do not bind in TF theory, because the errors in TF 
are far larger than the relatively small binding energies of molecules. 

6.4 Particles in boxes 

Before closing this chapter, let us return to the examples of the introduction and understand 
how they are connected to the Thomas-Fermi model. They highlight more clearly the relative 
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advantages and disadvantages of local approximations, since they apply in one dimension (no 
shell structure) and without interaction. 

First we note that, for the particles in a box, we cheated slightly, by applying the local 
approximation to the exact density. But if we're making the local approximation for the 
energy, where would we have gotten the exact density? A more honest calculation is to 
find the self-consistent density within the given approximation. Repeating the arguments of 
section 6.1, we find 

n(x) = 1 v^/i-VextO*)) (6.23) 

7T 

wherever the argument of the root is positive, zero otherwise, and p, is to be found by 
normalization (/ dxn(x) = N. 

For the particle in the box, v(x) vanishes inside the box, yielding 

N 

n (x) =n= — (6.24) 

Li 

This is quite different from the exact density for one electron, but the true density approaches 
it as N gets large. Using exactly this density, the energy contains only one term, the 
asymptotic one, E = 7r 2 iV 3 /(6L 2 ). 

Exercise 27 Energy from local approximation on exact density in box 

A 'fun ' math problem is to show that, using the local approximation evaluated on the exact 
density for a particle in a box, the energy can be found analytically to be 

You can check how it compares to your answers to Exercise X. 

In particular, for one electron, the self-consistent solution is too small by a factor of 3! 
Compare this with the 25% error using the exact density. This is because, although the local 
approximation is good for the energy, its derivative is not so good, and the density is quite 
sensitive to this. 

Note that our reasoning applies to all bounded Id problems, not just particles in boxes. 
Exercise 28 Large N for harmonic oscilator 

Calculate the self-consistent density of the local approximation for the harmonic oscillator, 
and compare with your exact densities. Also, compare the energies. 
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1. How would you go about finding the potential from the wavefunction of a two-electron 
system? For two electrons, can more than one potential have the same ground-state 
density? 

2. Define the ground-state wavefunction generating density n(r) without mentioning the 
external potential. 

3. When we say F[n] is a universal functional, exactly what do we mean? After all, since 
v ext (r) itself is a functional of n(r), is not the entire energy universal? 

4. Why do the HK theorems specify the ground-state and not, say, the first excited state? 

5. Using the local approximation for the one-dimensional kinetic energy developed in Chapter 
1, find the density for the one-dimensional harmonic oscillator, and evaluate the total 
energy on that density. Compare with your results for the HO problem there. 
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Kohn-Sham 



The most important step, establishing modern density functional theory as a useful tool, 
comes with the introduction of the Kohn-Sham equations, which allow almost all the 
kinetic energy to be calculated exactly. 

7.1 Kohn-Sham equations 

The major error in the Thomas-Fermi approach comes from approximating the kinetic energy 
as a density functional. We saw in chapter one that local approximations to the kinetic 
energy are typically good to within 10-20%. However, since for Coulombic systems, the 
kinetic energy equals the absolute value of the total energy, errors of 10% are huge. Even 
if we could reduce errors to 1%, they would still be too large. A major breakthrough in this 
area is provided by the Kohn-Sham construction of non-interacting electrons with the same 
density as the physical system, because solution of the Kohn-Sham equations produces the 
exact non-interacting kinetic energy, which includes almost all the true kinetic energy. 

We now have the theoretical tools to immediately write down these KS equations. Recall 
from the introduction, that the KS system is simply a fictitious system of non-interacting 
electrons, chosen to have the same density as the physical system. Then its orbitals are given 
by Eq. (1.5), i.e., 

{-^V 2 + ^(r)}^(r)= e ^;(r), (7.1) 



and yield 



N 



«(r) = EI*Wr (7-2) 



i=i 



The subscript s denotes single-electron equations. But the Euler equation that is equivalent 
to these equations is 

m +v - {I)=li > (7 - 3) 



where 



T s [n] = min($ T $), (7.4) 
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is the kinetic energy of non-interacting electrons. We have implicitly assumed that the Kohn- 
Sham wavefunction is a single Slater determinant, which is true most of the time. We denote 
the minimizing wavefunction in Eq. (7.4) by $[n]. This gives us a verbal definition of the 
Kohn-Sham wavefunction. The Kohn-Sham wavefunction of density n(r) is that wave- 
function that yields n(r) and has least kinetic energy. Obviously T s [n] = ($[n] T $[«-]), 
which differs from T[n]. 

Exercise 29 Kinetic energies 

Show that T[n] > T s [n}. What is the relation between T HF [n] and T s [n] ? Is there a simple 
relation between T[n] andT HF [n]? 

Now, write the ground-state functional of an interacting system in terms of the non- 
interacting kinetic energy: 

F[n] = T s [n] + U[n] + E xc [n]. (7.5) 

We have explicitly included the Hartree energy, as we know this will be a large part of the 
remainder, and we know it explictly as a density functional. The rest is called the exchange- 
correlation energy. Inserting F[n] into Eq. (6.13) and comparing with Eq. (7.3), we find 

Ws(r) = iw(r) + / d\^^- + v xc [n](r) v xc (r) = ^ (7.6) 
r — r ' dn(r) 

This is the first and most important relationship of exact density functional theory: from the 
functional dependence of F[n], we can extract the potential felt by non-interacting electrons 
of the same density. 

We note several important points: 

• The Kohn-Sham equations are exact, and yield the exact density. For every physical 
system, the Kohn-Sham alter ego is well-defined and unique (recall Fig. 1.2). There is 
nothing approximate about this. 

• The Kohn-Sham equations are a set of single-particle equations, and so are much easier 
to solve than the coupled Schrodinger equation, especially for large numbers of electrons. 
However, in return, the unknown exchange-correlation energy must be approximated. 
(We do not know this functional exactly, or else we would have solved all Coulomb- 
interacting electronic problems exactly.) 

• While the KS potential is unique by the Hohenberg-Kohn theorem (applied to non- 
interacting electrons), there are known examples where such a potential cannot be found. 
In common practice, this has never been a problem. 

• The great advantage of the KS equations over Thomas-Fermi theory is that almost all 
the kinetic energy (T s ) is treated exactly. 
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• The KS orbitals supercede the HF orbitals, in providing an exact molecular orbital theory. 
With exact KS theory, we see now how an orbital calculation can provide exact energetics. 
The HF orbitals are better thought of as approximations to exact KS orbitals. So the 
standard figure, Fig. 1.3, and other orbital pictures, are much more usefully interpreted 
as pictures of KS atomic and molecular orbitals. 

• Note that pictures like that of Fig. 1.2 tell us nothing about the functional dependence of 
u s (r). That KS potential was found simply by inverting the non-interacting Schrodinger 
equation for a single orbital. This gives us the exact KS potential for this system, but 
tells us nothing about how that potential would change with the density. 

• By subtracting T s and U from F, what's left (exchange and correlation) will turn out to 
be very amenable to local-type approximations. 

7.2 Exchange 

In density functional theory, we define the exchange energy as 

E x [n] = {<S>[n]\v ee \<S>[n]) - U[n], (7.7) 

i.e., the electron-electron repulsion, evaluated on the Kohn-Sham wavefunction, yields a direct 
contribution (the Hartree piece) and an exchange contribution. In most cases, the Kohn-Sham 
wavefunction is simply a single Slater determinant of orbitals, so that E x is given by the Fock 
integral, Eq. (5.13) of the KS orbitals. These differ from the Hartree-Fock orbitals, in that 
they are the orbitals that yield a given density, but that are eigenstates of a single potential 
(not orbital-dependent). Note that Eq. (5.13) does not give us the exchange energy as an 
explicit functional of the density, but only as a functional of the orbitals. The total energy in 
a Hartree-Fock calculation is extremely close to T s + U + V ext + E x . 

The differences between HF exchange and KS-DFT exchange are subtle. They can be 
thought of as having two different sources. 

1. The KS exchange is defined for a given density, and so the exact exchange of a system 
is the exchange of the KS orbitals evaluated on the exact density. The HF exchange is 
evaluated on the HF orbitals for the system. 

2. To eliminate the density difference, we can compare KS £ x [n HF ] with that from HF. The 
remaining difference is due to the local potential for the KS orbitals. 

Exercise 30 HF versus DFT exchange 

Prove E HF < T s [n HF ] + U[n HF ] + V ext [n HF ] + E x [n HF ] for any problem. 
Exercise 31 Exchange energies 

Calculate the exchange energy for the hydrogen atom and for the three-dimensional harmonic 
oscillator. If I construct a spin singlet, with both electrons in an exponential orbital, what is 
the exchange energy then? 
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atom 


E T V ext V ee 


T s U E x 


T c U c E c 


He 
Be 
Ne 


-2.904 2.904 -6.753 0.946 
-14.667 14.667 -33.710 4.375 
-128.94 128.94 -311.12 53.24 


2.867 2.049 -1.025 
14.594 7.218 -2.674 
128.61 66.05 -12.09 


0.037 -0.079 -0.042 
0.073 -0.169 -0.096 
0.33 -0.72 -0.39 



Table 7.1: Energy components for first three noble gas atoms, found from the exact densities. 



7.3 Correlation 

In density functional theory, we then define the correlation energy as the remaining unknown 
piece of the energy: 

E c [n] = F[n] - T s [n] - U[n] - E x [n}. (7.8) 

We will show that it usually better to approximate exchange-correlation together as a single 
entity, rather than exchange and correlation separately. Inserting the definition of F above, 
we find that the correlation energy consists of two separate contributions: 

E c [n]=T c [n] + U c [n] (7.9) 

where T c is the kinetic contribution to the correlation energy, 

T [n] = T[n] - T s [n] (7.10) 

(or the correlation contribution to kinetic energy), and U c is the potential contribution to the 
correlation energy, 

U [n] = V ee [n] - U [n] - E x [n]. (7.11) 
For many systems, T c ~ — E c ~ —U c /2. 
Exercise 32 Correlation energy from Hamiltonian 

Define the DFT correlation energy in terms of the Hamiltonian evaluated on wavefunctions. 

Exercise 33 Correlation energy 

Prove E c < 0, and say which is bigger: E c or _E' rad . 

Finally, to get an idea of how large these energies are for real systems, in Table 7.3 we 
give accurate values for three noble gas atoms. These numbers were found as follows. First, 
a highly accurate solution was found of the full Schrodinger equation for each atom. From 
it, the total ground-state energy and its various components could be calculated. Also, the 
ground-state density was extracted from the wavefunction. Next, a search was made for the 
unique KS potential that corresponded to that density. This can be thought of as guessing 
the potential, solving for its density, and comparing with the exact density. Then changing the 
potential to shift the computed density toward the exact one, and repeating until converged. 
Once the KS potential and its orbitals are found, it is straightforward to evaluate its kinetic, 
Hartree, and exchange energies. The last set of columns were found by subtracting KS 
quantities from their exact counterparts. 

There are many points to note in this table: 
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• The exact kinetic energy is simply \E\, by the virial theorem for atoms. 

• The interaction with the nucleus, V ex t, is a little more than twice the kinetic energy, as 
it must also overcome the (positive) V ee . 

• The electron-electron repulsion is typically less than half the kinetic energy, but all three 
are on the same scale. Recall from Thomas-Fermi theory that V ext = —TTj'i = — 7V ee . 

• The KS kinetic energy T s is almost as large as the true kinetic energy, and they differ by 
less than 2%. 

• The Hartree energy is typically a small overestimate of Vee- 

• The exchange energy cancels a fraction of the Hartree energy, that fraction getting smaller 
as Z increases. 

• The kinetic contribution to the correlation energy, which is also the correlation contribu- 
tion to the kinetic energy, is on the same scale as — E . 

• The potential contribution to the correlation energy, which is also the correlation contri- 
bution to the potential is energy, is a little more than — 2T C . 

• The magnitude of all energy components grows with Z. 

We end with an exercise designed to make you think about the trends. 
Exercise 34 Energy components for atoms 

Use the data in Table 7.3 to answer the following: 

1. How close are the ratios ofT to V ee to V ex t to their Thomas-Fermi ideal values? What 
is the trend with Z? Comment. 

2. What is the ratio of\E x \ to U? How does it change with Z? Comment. 

3. Repeat above forT c versus T s . Comment. 

4. What is the ratio T c versus \ U C \ ? As a function of Z? Comment. 

7.4 Questions 

All questions are conceptual. 

1. Define the ground-state Kohn-Sham wavefunction generating density n(r). 

2. Consider two normalized orbitals in one dimension, 4>i(x) and 4>2{x), and the density 
n = \<t>i\ 2 + W 2 - What is the Kohn-Sham kinetic energy of that density? How does 
it change when we alter one of the orbitals? Repeat the question for the Kohn-Sham 
exchange energy. 
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The exchange energy of the orbitals is the same in Hartree-Fock and Kohn-Sham theory. 
What is the relation between the two for the Hartree-Fock density? 

Considering the exact relations for the asymptotic decay of the density of Coulombic 
systems, is there any significance to enoMO for the KS system? (See section 14.2). 



Chapter 8 



The local density approximation 

1 This chapter deals with the mother of all density functional approximations, namely 
the local approximation, for exchange and correlation. 



8.1 Local approximations 



To illustrate why a local approximation is the first thing to try for a functional, we return 
to the problem in Chapter 2 of the perimeter of a curve r{&). Suppose we live in a flatland 
in which it is very important to estimate the perimeters of curves, but we haven't enough 
calculus to figure out the exact formula. 

To study the problem, we first consider a smoother class of curves, 

r{6) = r (l + e cos(ra0)). (8.1) 

where n = 1,2,3... is chosen as an integer to preserve periodicity. Here e determines how 
large the deviation from a circle is, while n is a measure of how rapidly the radius changes 
with angle. To illustrate, we show in Fig. 8.1 the case where n = 20 and e = 0.1. 



n=20, eps=0.1 




circle 






Figure 8.1: Smooth parametrized curve, r = 
1 ©2000 by Kieron Burke. All rights reserved. 
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l+O.lcos(2O0). 
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Now, in flatland, we wish to approximate the perimeter with a local functional of r, so we 
write 

P loc = f d6 f(r) (8.2) 

To determine the function /, we note that the perimeter of a circle is 2-7rr. The only place 
our approximation can be exact is for the circle itself, since its r is independent of 6. Thus 
we should choose / = r if we want our functional to be correct whenever it can be. Thus 

P loc = d6 r{6) (8.3) 

is the local approximation to the perimeter. Our more advanced mathematicians might be 
able to prove simple rules, such as P loc < P, etc. Also, we can test this approximation on 
the squares we really care about, and find, for the unit square, P loc = 41n(l + \/2), a 12% 
underestimate, just like the kinetic energy approximation from Chapter 1. 



8.2 Local density approximation 

So now its simple to see how to construct local density approximations. 

Recall the example of Chapter 1, the non-interacting kinetic energy of same-spin electrons, 
as a warm-up exercise. We first write 

T loc [n] = r dxtJn(x)) (8.4) 

Next we deduced that t s , the kinetic energy density, must be proportional to n'' , from dimen- 
sional analysis, i.e., 

t,(n) = a s n 3 (8.5) 

Last, we deduced a 3 by ensuring T g loc is exact for the only system for which it can be exact, 
namely a uniform gas of electrons. We did this by looking at the leading contribution to the 
energy of the electrons as the number became large. This yielded a s = ir 2 /6. 

At this point, we introduce a useful concept, called the Fermi wavevector. This is the 
wavevector of the highest occupied orbital in our system. Since we have N electrons in the 
box, 

Ntt 

k F = —— = nix (ID polarized) (8.6) 

Then 

T s loc = ! X dxn(x)T s (n{x)), r s (ra) = k%/6 (8.7) 

J—oc 

Thus, for 3D Coulomb-interacting problems, we write 

E%[n] = J d?r f(n(v)) (8.8) 

where f(n) is some function of n. To determine this function, we look to the uniform electron 
gas. 
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8.3 Uniform electron gas 

The local approximation is exact for the special case of a uniform electronic system, i.e., 
one in which the electrons sit in an infinite region of space, with a uniform positive external 
potential, chosen to preserve overall charge neutrality. The kinetic and exchange energies 
of such a system are easily evaluated, since the Kohn-Sham wavefunctions are simply Slater 
determinants of plane waves. The correlation energy is extracted from accurate Monte Carlo 
calculations, combined with known exact limiting values. 

If we repeat the exercise above in three dimensions, we find that it is simpler to use 
plane-waves with periodic boundary conditions. The states are then ordered energetically by 
momentum, and the lowest occupied levels form a sphere in momentum-space. Its radius is 
the Fermi wavevector, and is given by 

k F = (3Tr 2 n)^ 3 , (8.9) 

where the different power and coefficient come from the different dimensionality. Note that we 
have now doubly-occupied each orbital, to account for spin. Dimensional analysis produces 
the same form for the kinetic energy as in ID, but matching to the uniform gas yields a 
different constant: 

T s TF [n] = ^- fcfr 4(r)n(r) = A s ( d 3 r n 5/3 (r) (3D unpolarized) (8.10) 

1(J J ■* 

where A s = (3/10) (3tt 2 ) 2 / 3 = 2.871. 

For the exchange energy of the uniform gas, we simply note that the Coulomb interaction 
has dimensions of inverse length. Thus we know that its energy density must be proportional 
to k F , leading to 

E™ A [n] = A x f d 3 r n 4/3 (r) (8.11) 



Evaluation of the Fock integral Eq. (5.13) for a Slater determinant of plane-wave orbitals 
yields the exchange energy per electron of a uniform gas as 

4 nif (") = ^ (8-12) 

so that A x = -(3/4)(3/tt) 1 /-'' = -0.738. 

Correlation is far more sophisticated, as it depends explictly on the physical ground-state 
wavefunction of the uniform gas. Another useful measure of the density is the Wigner-Seitz 
radius 

/ 3 \ 1/3 1.919 ,„ „„. 

r -=te) (8 - 13) 

which is the radius of a sphere around electron such that the volume of all the spheres matches 
the total density of electrons. Thus r s — > is the high-density limit, and r s — > oo is the low 
density limit. 
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We thus write 

£ c LDA = J d 3 rn(r) e™ if (r,,(r)) (8.14) 

where e" mf (r,) is the correlation energy per electron of the uniform gas. A simple way to 
think about correlation is simply as an enhancement over exchange: 

eZ ! (r s ) = F xc (r s )eT i! (r s ) (8.15) 

and this enhancement factor is plotted in Fig. 15.1. There are several important features to 
the curve: 

• At r s = (infinite density), exchange dominates over correlation, and F x = 1. 

• As r s — > (high density), there is a sharp dive toward 1. This is due to the long-ranged 
nature of the Coulomb repulsion in an infinite system, leading to 

e™ if (r s ) -> 0.0311 In r s - 0.047 + 0.009r s In r s - 0.017r s (r s -> 0) (8.16) 

The logarithmic divergence is not enough to make correlation as large as exchange in this 
limit, but we will see later that finite systems have finite correlation energy in this limit. 

• In the large r s limit (low-density) 

er ! (r s ) - - J + % - .... (r. - oo) (8.17) 

where the constants are d = 0.896 and d\ = 1.325. The constant d Q was first deduced 
by Wigner from the Wigner crystal for this system. Note that this means correlation 
becomes (almost) as large as exchange here, and so F™ lf (r, — > oo) = 1.896. Note also 
that the approach to the low-density limit is extremely slow. 



I 1 1 1 r 




1 2 3 4 5 6 

Figure 8.2: Exchange and correlation energies per particle for a uniform electron gas. 

The uniform gas exchange and correlation energies/particle are plotted in Fig 8.2, as a 
function of r s . The exchange is very simple, being proportional to l/r s . The sharp downturn 
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in the correlation as r s — > is due to the logarithmic term (but is not as singular as exchange, 
which is diverging as l/r s ). At large r s , the two curves become comparable, but r s must be 
larger than shown here. For typical valence electrons, r 3 is between 1 and 6, and correlation 
is much smaller than exchange. 

Over the years, this function has become very well-known, by combining limiting informa- 
tion like that above, with accurate quantum Monte Carlo data for the uniform gas. An early 
popular formula in solid-state physics is PZ81, while the parametrization of Vosko-Wilkes- 
Nusair (VWN) has been implemented in quantum chemical codes. Note that in the Gaussian 
codes, VWN refers to an older formula from the VWN paper, while VWN-V is the actual 
parametrization recommended by VWN. More recently, Perdew and Wang reparametrized the 
data. These accurate parametrizations differ only very slightly among themselves. 

An early but inaccurate extrapolation from the low-density gas, was given by Wigner: 

e c (r s ) = -a/(b + r s ) (Wigner) (8.18) 

where a = —.44 abd b = 7.8. This gives a ballpark number, but misses the logarithmic 
singularity as r s — ► 0. 

Unfortunately, we must postpone a general review of the performance of LDA until after 
the next chapter. This is because we have only developed density functional theory so far, 
and not accounted for spin effects. A slight generalization of LDA, called the local spin- 
density approximation, is what is actually used in practice. In fact, many (if not most) 
practical problems require this. For example, pure LDA performs badly whenever there is an 
odd number of electrons, since it makes no distinction between polarized and unpolarized 
densities. We emphasize this fact in the next few exercises. 

Exercise 35 Hartree energy for an exponential 

1. For an exponential density, n(r) = \JZ 3 /Tiexp(—2Zr), calculate (a) the Hartree poten- 
tial and (b) the Hartree energy. 

2. Find the exact exchange energy for the Hydrogen atom. 

3. Find the exact exchange energy for the Helium atom, using the exponential approximation 
of Chapter 4. 

Exercise 36 LDA exchange energy for 1 or 2 electrons 

1. Calculate the LDA exchange energy for the Hydrogen atom, and compare with exact 
answer. 

2. Find the LDA exchange energy for the Helium atom, using the exponential approximation 
of Chapter X, and compare with exact answer. 
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8.4 Questions 

1. Following the previous question, can you deduce the asymptotic form of the exchange- 
correlation potential, i.e., t; xc (r) as r — > oc? Would a local approximation capture this 
behavior? (See section 14.2). 



Chapter 9 



Spin 



In practice, people don't use Kohn-Sham density functional theory, they use Kohn-Sham 
spin-density functional theory. For many problems, separate treatment of the up and 
down spin densities yields much better results. 

9.1 Kohn-Sham equations 

We next introduce spin density functional theory, a simple generalization of density functional 
theory. We consider the up and down densities as separate variables, defined as 

n a (r) = N J dx2 ■ ■ ■ J dxN |\P(r, a, X2, ■ ■ ■ , xn)\ 2 , (9-1) 

with the interpretation that n a (y)d?r is the probability for finding an electron of spin a 
in d 3 r around r. Then the Hohenberg-Kohn theorems can be proved, showing a one-to- 
one correspondence between spin densities and spin-dependent external potentials, v ext)lT (r). 
1 Similarly, the Kohn-Sham equations can be developed with spin-dependent Kohn-Sham 
potentials. All modern density functional calculations are in fact spin-density functional 
calculations. This has several advantages: 

1. Systems in collinear magnetic fields are included. 

2. Even when the external potential is not spin-dependent, it allows access to magnetic 
response properties of a system. 

3. Even if not interested in magnetism, we will see that the increased freedom in spin DFT 
leads to more accurate functional approximations for systems that are spin-polarized, 
e.g., the Li atom, whose spin densities are plotted in Fig. 9.1. 

We will see in the next sections that it is straightforward to turn some density functionals 
into spin-density functionals, while others are more complicated. 

1 This is not quite true. There is a little more wriggle room than in the DFT case. For example, for a spin-up H atom, the spin down potential 
is undefined, so long as it does not produce an orbital whose energy is below -1/2. 
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Figure 9.1: Radial spin densities in the Li atom. 



9.2 Spin scaling 



We can easily deduce the spin scaling of orbital functionals, while little is known of the exact 
spin-scaling of correlation in general. By orbital functionals, we mean those that are a sum 
of contributions from the orbitals of each spin seperately. 

We do spin-scaling of the kinetic energy functional as an example. Suppose we know T s [n] 
as a density functional for spin-unpolarized systems. Let's denote it T s unpol [n] to be very 
explicit. But since the kinetic energy is the sum of contributions from the two spin channels: 



T s [n h ni] = r s[«T>°] +T s [Q,ni] 



(9.2) 



i.e., the contributions come from each spin separately. Applying this to the spin-unpolarized 
case, we find 

T s ! mpol W = T s [n/2, n/2] = 2%\n/2, 0] (9.3) 
or T s [n, 0] = T s unpol [2n]/2. Inserting this result back into Eq. (9.2), we find 



T s [n hni ] = - {r & unpol [2n T ] + T^ 1 ^}} 



(9.4) 



This clearly yields a consistent answer for an unpolarized system, and gives the result for a 
fully-polarized system: 

T pol [n] = T s ,mpol [2n]/2 (9.5) 



Analogous formulas apply to exchange: 

E x [n hni ] = \ {£r ol [2n T J + E^n,]} 



and 



El°\n] = E™ pol {2n]/2 



.6) 
'•7) 
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9.3 LSD 

The local spin density (LSD) approximation is the spin-scaled generalization of LDA. As input, 
we need the energies per particle of a spin-polarized uniform gas. The formula for the exchange 
energy is straightforward, but the correlation energy, as a function of spin-polarization, for 
the uniform gas, must be extracted from QMC and known limits. 

Begin with exchange. In terms of the energy density of a uniform gas, Eq. (9.6) becomes: 

e x (n h ni ) = \ { e r ol (2n T ) + er ol (2n|)} (9.8) 
Now we know e™ lpol (n) = A x n i/3 , so 

e x (n T , n { ) = 2 1 / 3 A x [nf + nf 3 } (9.9) 
Now we introduce a new, useful concept, that of relative spin polarization. Define 

C = n l^ (9.10) 
n 

as the ratio of the spin difference density to the total density. Thus £ = for an unpolarized 
system, and ±1 for a fully polarized system. In terms of C, 

n T = (l + CW2, th = (1-CW2, (9.11) 
Inserting these definition into Eq. (9.9), we find, for the energy density: 

e x (n, = eT*°\n) t 1 + 4/3 + (1 ~ C) 4/a (9 . 12) 

Thus our LSD exchange energy formula is 

^..J - .4, J „W'» W WglfcCg (9 , 3) 

where £(r) is the local relative spin polarization: 

CM = " T(r) :; i(r) (9-14) 
n(r) 

Unfortunately, the case for correlation is much more complicated, and we merely say that 
there exists well-known parametrizations of the uniform gas correlation energy as a function of 
spin polarization, e" mf (n, £). In Fig. 9.2, we show what happens. When the gas is polarized, 
exchange gets stronger, because electrons exchange only with like spins. On the other hand, 
exchange keeps the electrons further apart, reducing the correlation energy. This effect is 
typical of most systems, but this has not been proven generally. 

Exercise 37 Relative spin polarization 
Sketch £(r) for the Li atom, as shown in Fig. 9.1. 

Exercise 38 Spin polarized Thomas-Fermi 

Convert the Thomas-Fermi local approximation for the kinetic energy in 3D to a spin-density 
functional. Test it on the hydrogen atom, and compare to the regular TF. 
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Figure 9.2: Effect of (full) spin polarization on exchange and correlation energies of the uniform electron gas. 

9.4 Questions 
All are conceptual. 

1. Why does LSD yield a more accurate energy for the Hydrogen atom than LDA? 

2. State in words the relation between v s [n](r) and v s „[n^, nj(r). When do they coincide? 

3. Deduce the formula of T s loc for unpolarized electrons in a one-dimensional box. 




Chapter 10 
Properties 



In this chapter, we survey some of the more basic properties that are presently being 
calculated with electronic structure methods, and include how well or badly LDA does for 
these. The local density approximation to exchange-correlation was introduced by Kohn 
and Sham in 1965, and (the spin- dependent generalization) has been one of the most 
succesful approximations ever. Until the early 90' s, it was the standard approach for all 
density functional calculations, which were often called ab initio in solid state physics. It 
remains perhaps the most reliable approximation we have. 



10.1 Total energies 

For atoms and molecules, the total exchange energy is typically underestimated by about 
10%. On the other hand, the correlation energy is overestimated by about a factor of 2 or 3. 
Since for many systems of physical and chemical interest, exchange is about 4 times bigger 
than correlation, the overestimate of correlation compliments the underestimate of exchange, 
and the net exchange-correlation energy is typically underestimated by about 7%. 

Exercise 39 Wigner approximation for He atom correlation energy: 
By approximating the density of the He atom as a simple exponential with Z e ff = 7/4 as we 
found in Chapter 4, calculate the IDA correlation energy using the Wigner approximation, 
Eq. (8.18). 





£ x 


E c 


E xc 


atom 


LDA exact error % 


LDA exact error % 


LDA exact error % 


He 
Be 
Ne 


-0.883 -1.025 -14 % 
-2.321 -2.674 -13 % 
-11.021 -12.085 -9 % 


-0.112 -0.042 167 % 
-0.225 -0.096 134 % 
-0.742 -0.393 89 % 


-0.996 -1.067 -7 % 
-2.546 -2.770 -8 % 
-11.763 -12.478 -6 % 



Table 10.1: Accuracy of LDA for exchange, correlation, and XC for noble gas atoms. 
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10.2 Densities and potentials 



The self-consistent LSD density is always close to the exact density, being hard to distinguish 
by eye, as seen in Fig. 10.1. However, when we calculate n(r), as defined earlier, we can 
12 




Figure 10.1: Radial density of the Ne atom, both exactly and from an LDA calculation 

see the difference in the asymptotic region: Note that the density integrates to 9.9 at about 
-1 i 1 1 1 1 




Figure 10.2: Logarithmic derivative of density for Ne atom, both exactly and in LDA. The curves are indistinguishable close 
to the nucleus, where they both reach -10. 

r = 2.68. So for all the occupied region, even n(r) is extremely accurate in LDA, and all 
integrals over the density are very accurate. For example, (1/r) is 31.1 exactly, but 31.0 in 
LDA. 

The Kohn-Sham potentials reflect this similarity. In Fig. 10.3, we show both the LDA 
and exact potentials, and how close they appear. However, note that they are dominated 
by terms that are (essentially) identical in the two cases. The largest term is the external 
potential, which is — 10/r. This determines the scale of the figure. The next term is the 
Hartree potential, plotted in Fig. 10.4. This is essentially the same in both cases, since the 
densities are almost identical. Again note the scale. 
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Figure 10.3: Kohn-Sham potential of Ne atom, both exactly and in LDA. 
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r 

Figure 10.4: Hartree potential of Ne atom, both exactly and in LDA. 

The LDA exchange-correlation potential, however, looks very different from the exact 
quantities for any finite systems, as shown in Fig. 10.5. This in turn means that the 
orbital eigenvalues can be very different from exact Kohn-Sham eigenvalues. Thus ionization 
potentials from orbital energy differences are very poor. This will be discussed in great detail 
in chapter X. 



10.3 Ionization energies and electron affinities 

In Table ??, we list ionization potentials for the first twenty atoms. There are many things 
to see in this table: 

• Comparing the exchange results (which are essentially identical to HF) with the exact 
results, we find that HF underestimates ionization potentials by about 1 eV . In fact, the 
mean error is -0.9 eV. 

• Repeating with the LSD numbers, we find that the errors vary in sign. Now we must 
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atom 


LSDX 


X 


LSD 


exact 


h 


12.40 


13.61 


13.00 


13.61 


he 


22.01 


23.45 


24.28 


24.59 


li 


5.02 


5.34 


5.46 


5.39 


be 


7.63 


8.04 


9.02 


9.32 


b 


7.52 


7.93 


8.57 


8.30 


c 


10.72 


10.81 


11.73 


11.26 


n 


13.91 


14.00 


14.93 


14.55 





11.78 


11.87 


13.85 


13.62 


f 


16.14 


15.66 


17.98 


17.45 


ne 


20.35 


19.81 


22.07 


21.62 


na 


4.84 


4.94 


5.32 


5.13 


mg 


6.48 


6.59 


7.70 


7.64 


al 


5.14 


5.49 


5.99 


5.99 


si 


7.40 


7.65 


8.26 


8.16 


P 


9.63 


10.04 


10.51 


10.53 


s 


8.79 


9.03 


10.52 


10.37 


cl 


11.67 


11.77 


13.22 


12.98 


ar 


14.42 


14.76 


15.90 


15.84 



Table 10.2: Ionization energies of first twenty atoms, non-relativistic, using exact-exchange densities, in eV. The X results 
use Engel's code, and the LSD results are evaluated on those densities. The 'exact' results are from Davidson. 
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r 

Figure 10.5: Exchange-correlation potential of Ne atom, both exactly and in LDA. 

take the mean absolute error (MAE), which is 0.25 eV. This is much better than HF, 
but not accurate enough for many purposes. 

• Interestingly, if we compare LSDX with the exact exchange numbers, we find that LSDX 
almost always underestimates the ionization energy, and has an MAE of 0.4 eV. Thus 
LSDX is a poorer approximation to exact exchange than LSD is to the exact result. 
Again, there is a cancellation of errors between exchange and correlation. 

Electron affinities 
Electronegativity and Hardness 

10.4 Dissociation energies 

We consider first atomization energies, being the energy difference between molecules or 
solids at equilibrium, and their constituent atoms. These are called cohesive energies of 
solids. To calculate these, the zero-point energy of the molecule (or solid) must be added to 
the minimum in its total energy curve. Typically, HF underestimates such binding energies 
by about 100 kcal/mol, while LDA overbinds by about 30 kcal/mol, or 1.3 eV, or about 50 
millihartrees. This meant LDA was never adopted as a general tool in quantum chemistry. For 
similar reasons, this also leads to transition state barriers that are too low, often non-existent. 

Finally, qualitative errors occur for highly correlated systems, such as the solid NiO or the 
molecules Cr 2 , or H 2 stretched to large distances. These systems are extremely difficult to 
study with density functional approximations, as will be discussed below. 

10.5 Geometries and vibrations 

However, bond lengths are extremely good in LDA, usually (but not always) being underes- 
timated by about 1-2%. 
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10.6 Transition metals 

In solid state physics, an infamous failure is the prediction that the non-magnetic structure 
of iron is of slightly lower energy than the magnetic one. 

In general, density functionals are much less reliable for transition metal complexes than 
first- and second-row elements, with bond-length errors more likely in the 0.05Arange. 

10.7 Weak bonds 

The remarks above apply to covalent, metallic, and ionic bonds. But many processes are 
determined by much weaker bonds, such as hydrogen bonds and dispersion forces. LDA 
significantly overbinds both hydrogen bonds, and van der Waals dimers. 

Between two isolated pieces of matter (no significant density overlap) there are always van 
der Waals forces due to fluctuating dipoles. In absence of permanent dipoles, and ignoring 
relativistic effects, the energy between two such pieces behaves as 

C 

E bi „d = --^, R^oo (10.1) 

where R is their separation. The exchange energy is always repulsive, so the existence of Ce 
is a purely correlation effect. 

LDA cannot reproduce this asymptotic behavior, since any contribution to the energy 
difference must come from a density difference. Thus the LDA correlation energy vanishes 
exponentially in this limit. 

10.8 Gaps 

On the other hand, for solids, these eigenvalues are often plotted as the band structure in 
solid-state texts. In these cases, the overall shape and position is good, but the band gap, 
between HOMO and LUMO, is consistently underestimated by at least a factor of 2. In some 
cases, some semiconductors have no gap in LDA, so it makes the incorrect prediction that 
they are metals. 

10.9 Questions 

1. Does the reliability of LSD show that most systems are close to uniform? 



Part III 
Analysis 
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Simple exact conditions 

1 In this chapter, we introduce various simple exact conditions, and see how well LDA 
meets these conditions. 

11.1 Size consistency 

An important exact property of any electronic structure method is size consistency. Imagine 
you have a system which consists of two extremely separated pieces of matter, called A and 
B. Their densities are given, n^fr) and n B (r), and the overlap is negligible. Then a size- 
consistent treatment should yield the same total energy, whether the two pieces are treated 
separately or as a whole, i.e., 

E[n A + n B ] = E[n A ]+E[n B }. (11.1) 

Configuration interaction calculations with a finite order of excitations are not size-consistent, 
but coupled-cluster calculations are. Most density functionals are size-consistent, but some 
with explicit dependence on N are not. Any local or semilocal functional is size-consistent. 

11.2 One and two electrons 

For any one electron system 

Ke = E x = -U E c = (N=l) (11.2) 

This is one of the most difficult properties for local and semilocal density functionals to get 
right, because the exact Hartree energy is a non-local functional. The error made for one- 
electron systems is called the self-interaction error, as it can be thought of as arising from 
the interaction of the charge density with itself in U. By spin-scaling, we find 

E x = -U/2 (N = 2, unpolarized) (11.3) 

for two spin-unpolarized electrons. LSD does have a self-interaction error, as can be seen for 
the H atom results in the tables. 

1 ©2000 by Kieron Burke. All rights reserved. 
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Exercise 40 Fermi- Amaldi correction 

An approximation for the exchange energy is E x [n] = —U[n]/N. This is exact for one or 
two (spin-unpolarized) electrons. Show that it is not size-consistent. 

This problem becomes acute when there is more than one center in the external potential. 
Consider, e.g. H \. As the nuclear separation grows from zero, the charge density spreads out, 
so that a local approximation produces a smaller and smaller fraction of the correct result. 

Exercise 41 Stretched 

Calculate the LSD error in the exchange energy of in the limit of infinite separation. 

11.3 Lieb-Oxford bound 

Another exact inequality for the potential energy contribution to the exchange-correlation 
functional is the Lieb-Oxford bound: 

U xc [n] > C L0 E^ DA [n} (11.4) 

where Clo < 2.273. It takes a lot of work to prove this result, but it provides a strong bound 
on how negative the XC energy can become. 

Exercise 42 Lieb-Oxford bound 

What does the Lieb-Oxford bound tell you about e£ mf (r s ) ? Does LDA respect the Lieb-Oxford 
bound? 

11.4 Bond breaking 

Another infamous difficulty of DFT goes all the way back to the first attempts to understand 
H 2 using quantum mechanics. In the absence of a magnetic field, the ground state of two 
electrons is always a singlet. At the chemical bond length, a Hartree-Fock single Slater deter- 
minant is a reasonable approximation to the true wavefunction, and yields roughly accurate 
numbers: 

$ HF (r l5 r 2 ) = 0(ri)0(r 2 ) « *(n, r 2 ) (11.5) 

However, when the bond is stretched to large distances, the correct wavefunction becomes 
a linear combination of two Slater determinants (the Heitler-London wavefunction), one for 
each electron on each H-atom: 

#(ri, r 2 ) = (0. 4 (ri)<2>B(r 2 ) + <M r i)<Mr 2 )) , (11.6) 

where 4>A is an atomic orbital centered on nucleus A and similarly for B. The total density is 
just n(r) = n^r) + n B (r), where n/(r) is the hydrogen atom density centered on nucleus I. 
But note that the density remains unpolarized. 
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If we apply LSD to the exact density, we make a large error, because E LSB [n] = E hDA [n A ] + 4 - what does the L bound say for one electron? 

E LDA [n B ], and we saw in a previous chapter the huge error LDA makes for the H-atom, 
because it does not account for spin-polarization. We would much rather find _E LSD = 

2E LSD (H). 

There is a way to do this, that has been well-known in quantum chemistry for years. If 
one allows a Hartree-Fock calculation to spontaneously break spin-symmetry, i.e., do not 
require the up-spin spatial orbital to be the same as that of the down-spin, then at a certain 
separation, called the Coulson-Fischer point, the unrestricted solution will have lower energy 
than the restricted one, and will yield the correct energy in the separated-atom limit, by 
producing two H atoms that have opposite spins. Thus we speak of RHF and UHF. The 
trouble with this is that the UHF solution is no longer a pure spin-eigenstate. Thus the 
symmetry dilemma is that an RHF calculation produces the correct spin symmetry with the 
wrong energy, while the UHF solution produces the right energy but with the wrong spin 
symmetry. 

Kohn-Sham DFT calculations with approximate density functionals have a similar problem. 
The LSD calculation will spontaneously break symmetry at a Coulson-Fischer point (further 
than that of HF) to yield the appropriate energy. Now the dilemma is no longer the spin 
eigenvalue of the wavefunction, but rather the spin-densities themselves. In the unrestricted 
LDA calculation, which correctly dissociates to two LSD H atoms, the spin-densities are 
completely polarized, which is utterly incorrect. 

11.5 Uniform limit 

We have already used the argument that any local approximation must use the value for the 
uniform limit as its argument. But more generally, we require 

E xc [n] -> Velf{n) n(r) -> n (constant) (11.7) 

where V is the volume of the system. This can be achieved, for example, by taking ever 
larger boxes, but keeping the density of particles fixed. This limit is obviously gotten right by 
LDA, but only if we use the correct energy density. Getting this right is relevant to accuracy 
for valence electrons of simple bulk metals, such as Li or Al. 

11.6 Questions 

1. If I approximate the correlation energy as about 1 eV per electron, is that a size-consistent 
approximation? 

2. A spin-restricted LSD calculation for two H atoms lOAapart yields a different answer to 
twice the answer for one H atom. Does this mean LSD is not size consistent? 

3. Does LDA satisfy the LO bound pointwise? Does LSD? 



Chapter 12 



Scaling 

1 In this chapter, we introduce the concept of scaling the density, as the natural way to 
understand the most important limits about density Junctionals. 

12.1 Wavefunctions 

It will prove very useful to figure out what happens to various quantities when the coordinates 
are scaled, i.e., when x is replaced by 72; everywhere in a problem, with 7 a positive number. 
For one electron in one dimension, we define the scaled wavefunction as 

7 (x) = 7^(71), < 7 < 00, (12.1) 

where the scaling factor out front has been chosen to preserve normalization: 

dxWx)\ 2 = dx^UHxM 2 = dx'Wx')\ 2 . (12.2) 

-00 J— 00 J— 00 

A scale factor of 7 > 1 will squeeze the wavefunction, and of 7 < 1 will stretch it out. 
For example, consider <p(x) = exp(— Then <f>~/(x) = ^/7exp(— -y |ar|) is a different 
wavefunction for every 7. For example, <t>2{x) = \/2exp(— 2|x|). In Fig. 12.1, we sketch 
a density for the one-dimensional H 2 molecule of a previous exercise, and the same density 
scaled by 7 = 2, 

n 7 (x) = 771(72:). (12.3) 
Note how 7 > 1 squeezes the density toward the origin, keeping the number of electrons 
fixed. We always use subscripts to denote a function which has been scaled. 
What is the kinetic energy of such a scaled wave-function? 

t[<m = ~£> iwi 2 = i£> i<^)i 2 = y £> >V)i 2 = i 2 m, (12.4) 

i.e., the kinetic energy grows quadratically with 7. The potential energy is not so straightfor- 
ward in general, because it is not a universal functional, i.e., it is different for every different 
problem. 

Vfa] = J^dxfaWVix) = f^dx'Mx'tfVi-) = (V(xh).) (12.5) 

1 ©2000 by Kieron Burke. All rights reserved. 
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Figure 12.1: Cartoon of the one-dimensional H 2 molecule density, and scaled by a factor of 2. 

However, if the potential energy is homogeneous of degree p, i.e., 

V(-yx) = l P V(x), (12.6) 

then 

V[cP 1 ]= 7 - p V[4>}. (12.7) 

Consider any set of trial wavefunctions. If scaling any member of the set leads to another 
member of the set, then we say that the set admits scaling. 

Exercise 43 Scaling of trial wavefunctions 

Do the sets of trial wavefunctions in Ex. 10 admit scaling? 

The only important change when going to three dimensions for our purposes in the pre- 
ceding discussion is that the scaling normalization factors change: 

<^ 7 (r) = 750(71-) n 7 (r) = 7 3 n(7r) (12.8) 

When including many electrons, one gets a similar factor for each coordinate: 

%( ri ,...,T N ) =7 3JV / 2 *( 7 r 1 ,...,7r J v) 1 (12.9) 

suppressing spin indices. We can they easily evaluate our favorite operators. It is simple to 
show: 

T[%] =7 2 r[^] (12.10) 

and 

v Be {%]=iv ee m (i2.il) 

while V ex t has no simple rule in general. 
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12.2 Density functionals 

Now we turn our attention to density functionals. Functionals defined as explicit operators 
on the Kohn-Sham wavefunction are far simpler than those including correlation effects. In 
particular, there is both the non-interacting kinetic energy and the exchange energy. We shall 
see that both their uniform scaling and their spin-scaling are straightforward. 

Consider uniform scaling of an iV-electron wavefunction, as in Eq. (12.9). The density of 
the scaled wavefunction is 

n(r) = N J d 3 r 2 ...J d 3 r N |* 7 (r, r 2 , . . . , r N )\ 2 , = 7 3 n( 7 r) = n 7 (r). (12.12) 

Now, a key question is this. If ^[n] is the ground-state wavefunction with density n(r), is 
^[n 7 ] = ^[n]? That is, is the scaled wavefunction the same as the wavefunction of the 
scaled density? We show below that the answer is no for the physical wavefunction, but is 
yes for the Kohn-Sham wavefunction. 

Consider the latter case first. We already know that T S [<J? 7 ] =7 2 T S [$]. Thus if $ minimizes 
T s and yields density n, then $ 7 also minimizes T s , but yields density n 7 . Therefore, $ 7 is 
the Kohn-Sham wavefunction for n 7 , or 

$ 7 [n] = $[n 7 ]. (12.13) 

This result is central to understanding the behaviour of the non-interacting kinetic and ex- 
change energies. We can immediately use it to see how they scale, since we can turn density 
functionals into orbital functionals, and vice versa: 

T s [n 7 ] = T s [$[n 7 ]] = T s [%[n}} = 7 2 T s [$[n]] = 7 2 T s [n], (12.14) 
using Eq. (12.24), and 

£ x [n 7 ] = £ x [$[n 7 ]] = E x [^[n]] = jE x [$[n}] = jE x [n], (12.15) 
using Eq. (F.15), and the fact that 

U[n 7 ] =jU[n] (12.16) 

These exact conditions are utterly elementary, saying only that if the length scale of a sys- 
tem changes, these functionals should change in a trivial way. Any approximation to these 
functionals must therefore satisfy these relations. However, they are also extremely powerful 
in limiting the possible forms functional approximations can have. 

The next exercise shows us that scaling relations determine the functional forms of the 
local approximation to these simple functionals, and the uniform electron gas is referred to 
only for the values of the coefficients. 

Exercise 44 Local density approximations for T s and E x : 

Show that, in making a local approximation for three-dimensional systems, the power of the 
density in both the non-interacting kinetic and exchange energy approximations are deter- 
mined by scaling considerations. 
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12.3 Correlation 

Levy has shown that for finite systems, E c [n n ] tends to a negative constant as 7 — » oc. From 
Fig. 12.2, we see that for the He atom E c varies little with scaling. We may write power 

He atom 



exact 

IDA ----- I 
1 23456789 10 

7 

Figure 12.2: Correlation energy of the He atom, both exactly and within LDA, as the density is squeezed. 

series for -E c [n 7 ] around the high-density limit: 

E c [n 7 ] = Ef ] [n] + E^[n]/"f + ..., 7 -> oc (12.17) 

We can easily scale any approximate functional, and so extract the separate contributions 
to the correlation energies, and test them against their exact counterparts, check limiting 
values, etc. A simple example is correlation in LDA. Then when n — ► n 7 , r s — > r s / 7 . Thus 

£c DA N = / d 3 r n 7 (r) e™«(r s ^r)) = f rf 3 rn(r) e r f (r,(r)/ 7 ) (12.18) 

In the high-density limit, 7 — > 00, e" mf ( r s/7) — * e c mf (0)- This is logarithmically singular in 
true LDA, as we saw in Chapter 8, and is an error made by LDA. 

From the figure, we see that for the He atom, as is the case with most systems of chemical 
interest, -E c [i"y] is almost constant. In fact, it is very close to linear in 1/7, showing that the 
system is close to the high-density limit. The high-density property is violated by LDA (see 
the figure), and diverges logarithmically as 7 — ► oc. This is another way to understand the 
LDA overestimate of correlation for finite systems. Since the LDA correlation energy diverges 
as 7 — > oc, it will be overestimated at 7 = 1. 

On the other hand, the correlation energy is believed to have the following expansion under 
scaling to the low-density limit: 

£ c [n 7 ] = 1 B[n] + j^ 2 C[n] + ..., 7^0 (12.19) 

Note that in this low-density regime, correlation is so strong it scales the same as exchange. 
This limit is satisfied within LDA, as seen from Eqs. (??), although it might not be numerically 
terribly accurate here. In the low density limit, e™ f (r s ) = —d /r s + d\jr'J 2 + . . ., so that 
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[n 7 ] = —jdol d 3 r n(r)/r s (r), vanishing linearly with 7, correctly. But most systems of 
chemical interest are closer to the high-density limit, where LDA fails. 

Exercise 45 Extracting exchange 

If someone gives you an exchange-correlation functional E xc [n], define a procedure for 
extracting the exchange contribution. 

As usual, one should not focus too much on the failings of the correlation energy alone. 
Remember, exchange scales linearly, and dominates correlation. When we add both contri- 
butions, we find Fig. 12.3. Clearly, the XC energy is best at 7 = 1, and worsens as 7 grows, 
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Figure 12.3: Exchange-correlation energy of the He atom, both exactly and within LDA. as the density is squeezed. 

because of the cancellation of errors. 

Exercise 46 Scaling LDA in Wigner approximation: 

Using the correct exchange formula and the simple Wigner approximation for the correlation 
energy of the uniform gas, Eq. (8.18), and the simple exponential density for the He atom, 
calculate the curves shown in Figs. 12.2 and 12.3. 

12.4 Correlation inequalities 

This reasoning breaks down when correlation is included, because this involves energies eval- 
uated on the physical wavefunction. As mentioned above, a key realization there is that the 
scaled ground-state wavefunction is not the ground-state wavefunction of the scaled density, 
because the physical wavefunction minimizes both T and V ee simultaneously. However, we 
can still use the variational principle to deduce an inequality. We write 

F[n,} = (*[n 7 ] | (f + V; e ) | *[ri,]) < [n] | (f + T/ ee ) | (12.20) 

or 

F[n-,} <7 2 7>]+7VeeM (12.21) 



98 CHAPTER 12. SCALING 

This is the fundamental inequality of uniform scaling, as it tells us inequalities about how 
correlation contributions scale. 

To see an example, apply Eq. (12.21) to n' = n 7 and write 7' = I/7, to yield 

F[n'] < T[n;,]/ 7 ' 2 + V„[n'^/i . (12.22) 

But we can just drop the primes in this equation, and multiply through by 7, and add 

(7 2 - 7)T[n] to both sides: 

7 2 T[n] + 7 Vee[n] < T[n 7 ]/ 7 + V ee {n.,\ + ( 7 2 - 7 )T[n]. (12.23) 

Now the left-hand-side equals the right-hand-side of Eq. (12.21), so we can combine the two 
equations, to find 

(7 - 1) T[ri n ] < 7 2 (7 - 1) T[n\. (12.24) 

Note that for 7 > 1, we can cancel 7 — 1 from both sides, to find T[n n ] < 7 2 T[n]. As 
we scale the system to high density, the physical kinetic energy grows less rapidly than the 
non-interacting kinetic energy. For 7 < 1, the reverse is true, i.e., 

T[n 7 ] < 7 2 T[n] 7 > 1, T[n 7 ] > 7 2 T[n] 7 < 1. (12.25) 

Exercise 47 Scaling V ee 

Show 

KeH > iVJn] 7 >1, KeK]<7KeN 7 < 1- (12.26) 

These inequalities actually provide very tight bounds on these large numbers. But of 
greater interest are the much smaller differences with Kohn-Sham values, i.e., the correlation 
contributions. So 

Exercise 48 Scaling E c 

Show 

E c [n 7 ] > jE c [n], 7 > 1, E c [n 7 ] < jE c [n], j < 1, (12.27) 

and 

T c [n 7 ] <j 2 T c [n], 7 > 1, T c [n 7 ] > j 2 T c [n], 7 < 1. (12.28) 
Exercise 49 Scaling inequalities for LDA 

Give a one-line argument for why LDA must satisfy the scaling inequalities for correlation 
energies. 

In fact, for most systems we study, both E c and T c are relatively insensitive to scaling the 
density toward the high-density limit, so these inequalities are less useful. 
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12.5 Virial theorem 

Any eigenstate wavefunction extremizes the expectation value of the Hamiltonian. Thus any 
small variations in such a wavefunction lead only to second order changes, i.e., 

E[4> s „i + Scj>] = E[<f> sol ] + 0(5<t> 2 ). (12.29) 

Assume <f> is the ground-state wavefunction of H. Then 7 is not, except at 7 = 1. But 
small variations in 7 near 1 lead only to second order changes in (H ), so 

d 

But, from Eq. (12.4), dT^/dj = 2jT{<j>], and from Eq. (12.5), dV[4)- f ]/dj = -j~ 2 (x % ) 
at 7 = 1, yielding 

2T=(x d ^). (12.31) 
dx 

The right-hand side is called the virial of the potential. So this must be satisfied by any 
solution to the Schrodinger equation, and can be used to test how accurately an approximate 
solution satisfies it. 

Note that this is true for any eigenstate and even for any approximate solution, once that 
solution is a variational extremum over a class of wave-functions which admits scaling, i.e., a 
class in which coordinate scaling of any member of the class yields another member of that 
class. Linear combinations of functions do not generally form a class which admits scaling. 

Exercise 50 Using the virial theorem 

Of the exercises we have done so far, which approximate solutions satisfy the virial theorem 
Eq. (12.31) exactly, and which do not? What can you deduce in the two different cases? 

We now derive equivalent formulae for the density functional case. We will first generalize 
Eq. (12.29), and then Eq. (12.31). As always, instead of considering the ground-state energy 
itself, we consider just the construction of the universal functional F[n]. Write F = T + V ee . 
Then 

F[n] = mn]\F\V[n]}, (12.32) 

and, if we insert any other wavefunction yielding the density n(r), we must get a higher 
number. Consider especially the wavefunctions * 7 [ni/ 7 ]. These yield the same density, but 
are not the minimizing wavefunctions. Then 

^| 7=1 <*> 4/7 ]|F|*> 1/7 ]) = 0. (12.33) 

But expectation values of operators scale very simply when the wavefunction is scaled, so 
that 

/l 7 =i [7 2 T[n lh ] + 7 Ke[m/ 7 ]) = 0, (12.34) 
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EAn]+TM = -^P 1 \ 1= i (12-36) 



or, swapping 7 to 1/7 in the derivative, 

2T[n] + V ee [n] = ^-| 7=1 (T[n 7 ] + KeM) . (12.35) 

This then is the generalization of Eq. (12.29) for density functional theory: Because ^>[n] 
minimizes F, variations induced by scaling (but keeping the density fixed) have zero derivative 
around the solution, leading to Eq. (12.35). 

12.6 Kinetic correlation energy 

As usual, the more interesting statement comes from the KS quantities. We apply Eq. (12.35) 
to both the interacting and Kohn-Sham systems, and subtract, to find: 

d E xc [n y ] t 
d"/ 

Exercise 51 Virial expressed as scaling derivative 
Prove Eq. (12.36). 

This statement is trivially true for Hartree and Exchange energies, so we subtract off the 
exchange contribution, and write: 

dE c [n 7 ] . 
d 7 

This is an extremely useful and powerful statement. It tells us how, if we're given a functional 
for the correlation energy, we can extract the kinetic contribution simply by scaling. 

Exercise 52 LDA kinetic correlation energy: 

Assuming the correlation energy density of a uniform gas is known, e" mf (r s ), give formulas 
for both t™ if (r s ) and < mif (r s ). 

Exercise 53 Scaling to find correlation energies: 
1. By applying Eq. (12.37) to n 7 (r), show 



T c [n] = ^P^| 7=1 - E c [n] (12.37) 



2. Next, consider Eq. (12.38) as a first-order differential equation in 7 for E [n~y\. Show 
that 

^T [ny], (12.39) 

using the fact that E c [nj] vanishes as 7 — > 0. 
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3. Show 

U c {n,]=2E c {n,]~ 1 dE j^ (12.40) 

and 

E c [n 1 ]= 1 2 p^U^n,,]. (12.41) 

4. Combine Eqs. (12.40) and (12.39) to find differential and integral relations between U c 
andT c . 

Exercise 54 Scaling of potentials: 

lfv c [n](r) = SE c [n]/n(r), and I define E^n]^) = E c [n 7 ], /)cw/su c [n](7, r) = 5E c [n](y)/n(v 
related to v c [n](r) ? 

Exercise 55 Scaling LDA in Wigner approximation: 

Within the simple Wigner approximation for the correlation energy of the uniform gas, 
Eq. (8.18), deduce both the potential and kinetic contributions, ?i" mf (r s ) and t" nlf (r s ), 
respectively. Then check all 6 relations of the previous exercise. 

12.7 Potential 

After our earlier warm-up exercises, it is trivial to derive the virial theorem for the exchange- 
correlation potential. For an arbitrary system, the virial theorem for the ground-state yields: 

2T= jr ($| ri . ViV"!*), (12.42) 

no matter what the external potential. For our problems, the electron-electron repulsion is 
homogeneous of order -1, and so 

2T + V ee = J d 3 r n(r)r ■ V« ex t.(r). (12.43) 

This theorem applies equally to the non-interacting and interacting systems, and by subtract- 
ing the difference, we find 

E xc [n] + T c [n] = - J d 3 r n(v) r ■ V?; xc (r). (12.44) 

Since we can either turn off the coupling constant or scale the density toward the high-density 
limit, this applies separately to the exchange contributions to both sides, and the remaining 
correlation contributions: 

E c [n] + T c [n] = - J d 3 r n(v) r ■ Vv c (r). (12.45) 

For almost any approximate functional, once a Kohn-Sham calculation has been cycled to 
self-consistency, Eq. (12.44) will be satisfied. It is thus a good test of the convergence of 
such calculations, but is rarely performed in practice. 
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Exercise 56 Virial for LDA: 

Show that the XC virial theorem is satisfied in an LDA calculation. 

12.8 Questions 

1. What is the relation between scaling and changing Z for the H-atom? 

2. Is the virial theorem satisfied for the ground-state of the problem V(x) = — exp(— |x|)? 

3. For a homogeneous potential of degree p, e.g., p=-l for the delta function well, we know 
E = —T. Can we find the ground-state by simply maximizing Tl 

4. What is the exact kinetic energy density functional for one electron in one-dimension? 
How does it scale? 

5. Does T s loc [n] from chapter ?? satisfy the correct scaling relation? (Recall question 1 from 
chapter ??). 

6. If W is the ground-state for an interacting electronic system of potential v ext (r), of what 
is * 7 a ground-state? 
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Adiabatic connection 



1 In this chapter, we introduce an apparently new formal device, the coupling constant of 
the interaction for DFT. This differs from that of regular many-body theory. But we see 
at the end that this is very simply related to scaling. 

13.1 One electron 

Introduce a parameter in H , say A. Then all eigenstates and eigenvalues depend on A, i.e., 
Eq. (1.9) becomes: 

H x d> x = eV- (13-1) 

Most commonly, 

H x = f + XV, 0<A<oc. (13.2) 

A A in front of the potential is often called the coupling constant, and then all eigenvalues 
and eigenfunctions depend on A. For example, for the harmonic oscillator of force constant 
k = uj 2 , the ground-state wavefunction is 4>(x) = (u;/7r) 1//4 exp(— cjx 2 /2). Since V(x) = 
kx 2 /2, \V(x) = Xkx 2 /2, i.e., the factor of A multiplies k. Then ui — ► \f\uj, and so 

fVx \ 1/4 

^(x) = ^ exp(-v / A^x 2 /2) (h armonic oscillator). (13.3) 

For a A-dependent Hamiltonian, we can differentiate the energy with respect to A. Be- 
cause the eigenstates are variational extrema, the contributions due to differentiating the 
wavefunctions with respect to A vanish, leaving 

Thus 

E = E x = l = E x =° + J Q d\{ X | | <p x ). (13.5) 



1 ©2000 by Kieron Burke. All rights reserved. 
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Furthermore, if A simply multiplies V, then 

j rpX 

^ r = (0 A |l>l^ A >, (13-6) 

so that 

E = E X=0 + f^d\{(j) X \V\(j) X ). (13.7) 
This is often called the Hellmann-Feynman theorem. 

Exercise 57 Hellmann-Feynman theorem 

Show that the 1-d H-atom and 1-d harmonic oscillator satisfy the Hellman-Feynman theorem. 
How does the approximate solution using a basis set do? 

In modern density functional theory, one of the most important relationships is between 
the coupling constant and coordinate scaling. We can see this for these simple 1-d problems. 
The Schrbdinger equation at coupling constant A is 

(f + XV(x)j <i> x {x) = E x <p x (x) (13.8) 

If we replace x by jx everywhere, we find 

(\f + A%z)) A ( 7 x) = E x <p x (jx) (13.9) 

Furthermore, if V is homogeneous of degree p, 

(f + \~/ 2 YV(x)) <t>\-/x) = 1 2 E X (f > x { 1 x) (13.10) 

Thus, if we choose \ ^i p+2 = 1, then Eq. (13.10) is just the normal Schrodinger equation, 
and we can identify (t> x (^x), which is equal (up to normalization) to <j> x {x) with <j>(x). If we 
scale both of these by I/7, we find 

(j) x (x) = 0i/ 7 (x) = A i/( P +2) (x) (13.11) 

In this case, changing coupling constant by A is equivalent to scaling by A 1 / 4 , but this 
relation depends on the details of the potential. 

Exercise 58 A-dependence 

Show that, for a homogeneous potential of degree p, 

E X = \^E (13.12) 
Exercise 59 A-dependence of 1-d H 

Find the coupling-constant dependence for the wavefunction and energy for the 1-d H-atom, 
V(x) = — 5(x), and for the 1-d harmonic oscillator, V(x) = \x 2 . 
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13.2 Adiabatic connection formula 

Now we apply the same thinking to DFT. When doing so, we always try to keep the density 
fixed (or altered only in some simple way), and not include the external potential. The 
adiabatic connection formula is a method for continuously connecting the Kohn-Sham system 
with the physical system. Very simply, we introduce a coupling constant A into the universal 
functional, to multiply V ee : 

F x \n] = min(*|f + A"l> ee I *) (13.13) 

For A = 1, we have the physical system. For A = 0, we have the Kohn-Sham system. We 
can even consider A — > oc, which is a highly correlated system, in which the kinetic energy is 
negligible. But for all values of A, the density remains that of the physical system. Note that 
this implies that the external potential is a function of A: f A x t( r )- We denote ^ x [n] as the 
minimizing wavefunction for a given A. We can generalize all our previous definitions, but we 
must do so carefully. By a superscript A, we mean the expecation value of an operator on 
the system with coupling-constant A. Thus 

T x [n] = (* A [n]|T|* A [n]), (13.14) 

but 

KeN = <* A M|AVy* A [n]>. (13.15) 
The Kohn-Sham quantities are independent of A, so 

E x c [n] = <* A [n] |f + Al/ ee | * A [n]) - <$[n] \f + XV ee \ <f[n}) (13.16) 

Exercise 60 Coupling constant dependence of exchange 
Show that 

U x [n] = \U[n], E x [n] = XE x [n], (13.17) 

i.e., both Hartree and exchange energies have a linear dependence on the coupling-constant. 

Our next step is to write the Hellmann-Feynman theorem for this A-dependence in the 
Hamiltonian. We write 

F[n] = T s [n] + jT 1 dX (9 x [n] |t> ee | * A [n]) (13.18) 

where F[n] = F x=1 [n], T s [n] = F x=0 [n}. The derivative w.r.t. A of F x [n] is just the 
derivative of the operator w.r.t. A, because $ A [n] is a minimizing wavefunction at each A, 
just as in the one-electron case. Inserting the definition of the correlation energy, we find 

Exc[n] = J Q [ dX <* A [n] |f ee | * A [n]) - 

= L id V U *M=J 1 dXU xc [n}W, (13-19) 
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where we have used = A(* A [«] |v^e| , 3 >A [n,]) and U x = XU, and introduced the simple 
notation U XC (X) = U xc /X for later convenience. This is the celebrated adiabatic connection 
formula. We have written the exchange-correlation energy as a solely potential contribution, 
but at the price of having to evaluate it at all intermediate coupling constants. We will see 
shortly that this price actually provides us with one of our most important tools for analyzing 
density functionals. 

A cartoon of the integrand in the adiabatic connection formula for a given physical system 
is sketched in Fig. 13.1. The curve itself is given by the solid line, and horizontal lines have 
been drawn at the value at A = and at A = 1. We can understand this figure by making 
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A 

Figure 13.1: Cartoon of the adiabatic connection integrand. 

several key observations: 

• The value at A = is just E x . 

• The value at A = 1 is just E x + U c . 

• The area between the curve and the x-axis is just E xc . 

• The area beneath the curve, between the curve and a horizontal line drawn through the 
value at A = 1, is just T c . 

Thus the adiabatic connection curve gives us a geometrical interpretation of many of the 
energies in density functional theory. 

To demonstrate the usefulness of adiabatic decomposition, we show this curve for the He 
atom in Fig. 13.2. Note first the solid, exact line. It is almost straight on this scale. This 
is telling us that this system is weakly correlated, as we discuss below. The numbers, from 
Table X, are E x = -1.025, E c = -0.042, U c = -0.079, and T c = 0.037. Thus E c is only 
slightly more negative than — T c . 

This produces yet another way to understand the cancellation of errors. We see that LDA 
improves as A grows from to 1. This is characteristic of a very general trend. But once we 
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Figure 13.2: Adiabatic decomposition of exchange-correlation energy in He atom, both exactly and in LDA. 

claim that the LDA error drops with A, then the usual cancellation of exchange and correlation 
errors follows, since the area under the curve will then be more accurate than the value at 
A = 0. Why this happens will be discussed later. But for now, we simply note that the 
statement that LDA improves with A is equivalent to the cancellation of errors statement, 
but is much more transparent. 

We can use the adiabatic decomposition (A-dependence of U XC (X)) to analyze any approx- 
imate functional, and find out why it behaves as it does. The next section shows how easy 
it is to find U XC (X) from E xc . In fact, we already know. 



13.3 Relation to scaling 

A third important concept is the relation between scaling and coupling constant. Consider 
<]> A [n], which minimizes T+XV ee and has density n(r). Then 9 x [n] minimizes T"/7 2 +(A/7)t4e 
and has density n 7 (r). If we choose 7 = 1/A, we find ^^[n-] minimizes A 2 (T + V ee ) and 
has density n 1 / A (r). But if A 2 (T + V^, e ) is minimized, then T + V ee is minimized, so that we 
can identify this wavefunction being simply 'I'Ini/.x]. If we scale both wavefunctions by A, we 
find the extremely simple but important result 

= *A[m/A], (13.20) 

which tells us how to construct a wavefunction of coupling constant A by first scaling the 
density by 1/A, finding the ground-state wavefunction for the scaled density, and then scaling 
that wavefunction back to the original size. 

This solves a mystery from the previous chapter: The ground-state wavefunction of a 
scaled density is the scaled ground-state wavefunction, but only if the coupling constant is 
changed, i.e., 

*[n 7 ] = *;/>]. (13.21) 
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Our argument also yields a simple relation for the energies: 

F x [n] = \ 2 F[n 1/x ], (13.22) 

which shows how the A-dependence of the universal functional (or, as we shall see, any energy 
component) is completely determined by its scaling dependence. Contrast Eq. (13.22) with 
Eq. (12.21). 

These relations prove to be extremely useful in analyzing functionals and their behavior. 
First note that we can take any functional and find out its scaling behavior quite easily. Then, 
through Eq. (13.22), applied to that functional, we can now deduce its coupling-constant 
dependence. For example, the adiabatic connection integrand of Eq. (13.19) is simply 

U xc {X)[n] =At/ xc [n 1/A ]. (13.23) 

The simplest example in this regard is the non-interacting kinetic energy, which we know 
scales quadraticly with scaling parameter. Then 

T s x [n] = X 2 T s [n i/x ] = T s [n] (13.24) 

i.e., the kinetic energy is independent of coupling constant, as it should be. Similarly, 

E x [n] = X 2 E x [n 1/x \ = XE x [n] (13.25) 

which again makes sense, since E x is constructed from the A-independent Kohn-Sham orbitals, 
integrated with A/|r — r'|. 

Much less trivial is the relation between different components of the correlation energy. 
Consider Eq. (13.19) in differential form: 

= <* A [n] \te\ *>]> - U[n] (13.26) 

Then, from Eq. (13.22), 

KeM = A(* A [n] |l/ ee | = X(U[n] + E x [n}) + U x [n], (13.27) 

where we have used the linear dependence of exchange on A, and we have written 

E c [n] = (T[n] - T s [n]) + (V ee [n] - U[n]) = T c [n] + U c [n] (13.28) 

and T c [n] is the kinetic contribution to the correlation energy, while U c [n] is the potential 
contribution. Inserting Eq. (13.27) into Eq. (13.26) and cancelling the trivial exchange 
contributions to both sides, we find 

= u x H/X (13.29) 

Finally, we rewrite this as a scaling relation. Since E x [n] = X 2 E c [ni/x] and U x [n] = 
X 2 U c [ni/x] and writing 7 = 1/A, we find, after a little manipulation 

07 
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which is just Eq. (12.38). Thus, the adiabatic connection formula may be thought of as an 
integration of this scaling formula between 7 = 1 and 7 = 00. 

Exercise 61 E c from T c 

Derive a relation to extract E c [n] from T c [n 7 ] alone. Rewrite this to get E*{n] from T£[n]. 

Exercise 62 Use Eq. (16.5) and the scaling relations of the previous chapter to write E c in 
terms of density matrices of different A. 

We can easily scale any approximate functional, and so extract the separate contributions 
to the correlation energies, and test them against their exact counterparts, check limiting 
values, etc. A simple example is LDA. Then when n — > n 7 , r s — > r s /j. Thus 

£ C LDA H = / rfV n 7 (r) er f (r s , 7 (r)) = / d 3 rn(r) e r f (r s (r)/ 7 ) (13.31) 

In the high-density limit, 7 — > 00, e" nlf (r s / 7 ) — ► e™ lf (0). This is singular in true LDA, and is 
an error made by LDA. However, in the low density limit, e" mf (r s ) = —d /r s + d[/r^ 2 + . . ., 
so that ££ DA [7i 7 ] = —jd J d A r n(r)/r s (r), vanishing linearly with 7, correctly. Note that in 
this low-density regime, correlation is so strong it scales the same as exchange. 

Exercise 63 Adiabatic connection for He atom: 

Using the Wigner approximation and the effective exponential density, draw the adiabatic 
connection curve for He. 

Exercise 64 Adiabatic connection for H atom: 



Draw the adiabatic connection curve for the H atom. 

Exercise 65 Changing A: 

Derive E^ A ' X \n] and T^ k \n]. 



13.4 Static correlation 

We have seen how, for the total XC energy of many systems of interest, the adiabatic 
connection curve is close to linear. In fact, it has always been found that the adiabatic 
connection curve is (slightly) concave upwards: 

d > {unproven) (13.32) 

d\ 

Given the geometric interpretation of Fig X, this means T c < \E C \. To quantify the fraction 
of correlation that is kinetic, we define 
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Thus, it has always been found that < b < 1/2, but b is always close to 0.5 for real systems. 
An alternative way to write it is 

E c = (1 - b)U c (13.34) 

or that if we simplify the adiabatic curve as a step down at some value of A from E x to 
E x + U c , then the step occurs at A = b. 

Considering both the figures and the numbers, we see that LDA clearly significantly un- 
derestimates b. This is again because of the logarithmic divergence as r s — ► 0, and we will 
see that this does not happen for more sophisticated approximations. 

The reason that most systems have b close to 1/2 is that they are not very far from 
the high-density limit. As one moves towards lower densities (much lower than in realistic 
systems), b reduces, even coming close to zero. Thus, even for the uniform gas, T c eventually 
becomes small relative to \U C \, as r s —> 00. Thus we can speak of b as measuring the amount 
of dynamic correlation, i.e., the fraction of correlation energy that is kinetic, with a value 
of 1/2 being the maximum. Typical systems are close to this value, and only for very low 
densities does b become small. 
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Figure 13.3: Adiabatic decomposition of exchange-correlation energy in stretched H2. 

To appreciate why we are introducing this parameter, we next consider what happens to 
a H 2 molecule when we stretch it. In Figure 13.3, we plot the adiabatic connection curve 
at R = 5. It looks very different from that of He, or of H2 at equilibrium. To understand 
it, we note that at A = 1, the curve is almost flat, and just about equal to —5/8, i.e., the 
exchange(-correlation) energy of two separate H atoms. On the other hand, at A = 0, we 
have the exact exchange value of KS theory. For this problem, the electrons are always in a 
singlet and the exchange energy is about -0.42. The value of b is about 0.16, and this system 
has strong static correlation. 

There is no way for LDA to produce a curve that looks anything like this. Because the 
average density will be almost the same as a single H atom, the adiabatic curve is very flat, 
with essentially no static correlation. Since we are beyond the Coulson-Fischer point, we 
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must chose if we wish to break spin-symmetry or not. If we allow symmetry breaking, we will 
get the A = 1 end about right, but have a flat curve, and so overestimate the magnitude of 
the correlation energy. On the other hand, if we don't, we will get a curve that bends just a 
little down from the A = end, and so underestimate the correlation. 

As we will see throughout the rest of this book, this is a very important problem for DFT 
in quantum chemistry. Many molecules exhibit similar behavior, not in their total energies, 
but in their energy differences between bonded and dissociated states. 

Exercise 66 Adiabatic connection for stretched H 2 

Draw an accurate adiabatic connection curve for H 2 with a bond length of 20 A. 

13.5 Questions 

1. Using the high-density limit of Wigner correlation, estimate roughly how large the corre- 
lation energy usually is within LDA. 
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Chapter 14 



Discontinuities 



14.1 Koopman's theorem 

Since the large distance condition Eq. (5.18) is true for both the interacting wavefunction 
and an independent particle description, the following approximate theorem was first noted 
by Koopmans (who later won a Nobel prize in economics). In a Hartree-Fock calculation, if 
you ignore relaxation, meaning the change in the self-consistent potential, when an electron 
is removed from the system, you find 

7 HF = E RF (N - 1) - E nF (N) « -e N (14.1) 

where en is the eigenvalue of the highest occupied orbital. 

Exercise 67 Koopman's theorem for Id He 
For the accurate HF 1-d He calculation, test Koopman's theorem. 

But, in exact density functional theory, Koopman's theorem is exact. To see why, note 
that Eq. (5.18) is true independent of the strength of the interaction. Consider the coupling 
constant A, keeping the density fixed. Since the density is fixed, its asympototic decay is the 
same for any value of A, i.e., the ionization potential is independent of A. In particular, at 
A = 0, its just -cHOMO^ w h ere we now re f er t0 tne HOMO of the KS potential. Thus 

/ = E(N - 1) - E(N) = - e H0M0 . (14.2) 

We will see below that this condition is violated by all our approximate functionals. 

Exercise 68 Asymptotic behavior of the density 

From the information in this and earlier chapters, deduce the limiting values of the asymptotes 
on the right of Fig. 5.2. 

14.2 Potentials 

It is straightforward to argue for the exact asymptotic behaviour of the exchange-correlation 
potential for a Coulombic system. Consider first exchange. Far from the nucleus, and 
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electron in an orbital should see an effective nuclear charge of Z — (N — 1), the TV — 1 
because there are N — 1 electrons close to the nucleus. Since v ext (r) = —Z/r, and at 
large distances v H (r) = Njr, this implies v x (r) = —l/r. Since this effect will occur in 
any wavefunction producing the exact density, this must also be true for v xc (r), so that the 
correlation contribution must decay more rapidly. In fact, 

v x (r) -> -l/r, v c (r) -> -a(N - l)/2r 4 (r -> oc) (14.3) 

where a(N — 1) is the polarizability of the (TV — l)-electron species. The decay of the 
correlation potential is so rapid that it has only ever been clearly identified in an exact 
calculation for H . 

Exercise 69 Asymptotic potentials 

Derive the asymptotic condition on v x (r) from the exact decay of the density. How does 
v x DA (r) behave at large distances? 

As mentioned earlier, LDA does not provide very realistic looking exchange-correlation 
potentials. Because the density decays exponentially in the tail region, so too will the LDA 
potential. In Fig. 14.1, we illustrate this effect. Of course, LDA also misses the shell 
structure between the Is and 2s and 2p electrons, just as smooth approximations for T s miss 
shell structure. This is partially explained by the above argument: the long-range decay of 




r 

Figure 14.1: XC potential in the Ne atom, both exact and in LDA 

v x (r) is due to non-self-interaction in the exact theory, an effect missed by LDA. This has a 
strong effect on the HOMO orbital energy (and all those above it), as will be discussed in 
more detail below. 

14.3 Derivative discontinuities 

We can understand the apparent differences between LDA and exact potentials in much more 
detail. Consider what happens when two distincct subsystems are brought a large distance 
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from each other. For example, take the H-He + case. At large distances, the interacting 
wavefu notion for the two electrons should become two (unbalanced) Slater determinants, 
making up a Heitler-London type wavefunction, yielding a density that is the sum of the 
isolated atom densities, centered on each nucleus. In our world of 1-d illustrations, this will 
be 

n(x) = exp(-2|.x|)+2exp(-4|x-L|) (14.4) 

where the H-atom is at and the He ion is at L. For this two electron system, we find 
the molecular orbital as (p(x) = \Jn(x)/2 and the Kohn-Sham potential by inversion. This 
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"Id H-He+ potential 



Figure 14.2: Density and Kohn-Sham potential for Id H-He + . 

is shown in Fig. 14.2, for L = 4. The apparent step in the potential between the atoms 
occurs where the dominant exponential decay changes. This is needed to ensure KS system 
produces two separate densities with separate decay constants. Elementary math then tells 
you that the step will be the size of the difference in ionization potentials between the two 
systems. 

More generally, for e.g., NaCI being separated, the step must be the difference in l-E.A., 
where E.A. is the electron affinity. This is the energy cost of transferring an infinitesimal of 
charge from one system to another, and the KS potential must be constructed so that the 
molecule dissociates into neutral atoms. But the same reasoning can be applied to a single 



116 



CHAPTER 14. DISCONTINUITIES 



atom, in contact with a bath of electrons, eg a metal. One finds that the step in the potential 
induced by adding an infinitesimal of charge is l-EA. 

But LDA is a smooth functional, with a smooth derivative. Thus its potential cannot 
discontinuously change when an infinitesimal of charge is added to the atom. So it will be in 
error about (l-EA)/2 for the neutral, and the reverse for the slightly negatively charged ion. 
Comparing the potentials for the Ne atom, we see in Fig. 14.3 how close the LDA potential 
is in the region around r = 1, relevant to the 2s and 2p electrons. However, it has been 
shifted down by -0.3, the LDA error in the HOMO orbital energy. Since Ne has zero electron 
affinity, this is close to (l-EA)/2 for this system, since 1=0.8. This rationalizes one feature of 
the poor-looking potentials within LDA. 




r 

Figure 14.3: XC potential in the Ne atom, both exact and in LDA, in the valence region, but with LDA shifted down by 
-0.3. 

If the functional is generalized to include fractional particle numbers via ensemble DFT, 
these steps are due to discontinuities in the slope of the energy with respect to particle number 
at integer numbers of electons. Hence the name. 

14.4 Questions 

1. What kind of chemical systems will suffer from strong self-interaction error? 

2. Which is better, to get the right energy and wrong symmetry, or vice versa? 

3. Is Koopmans' theorem satisfied by LDA? 

4. If one continues beyond x = 5 in Fig. 14.2, what happens? 
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Analysis tools 



In this chapter, we introduce a variety of tools that help analyze how approximate density 
junctionals work. 



15.1 Enhancement factor 
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Figure 15.1: Enhancement factor for correlation in a uniform electron gas as a function of Wigner-Seitz radius for unpolarized 
(solid line) and fully polarized (dashed line) cases. 



15.2 Density analysis 

What we can do is ask how accurately we need to know the uniform gas inputs. Recall the 
radial density plot of the Ar atom, Fig. 5.1. Now we plot the local Wigner-Seitz radius, 
r s (r), in Fig. 15.2, and find that the shells are less obvious. The core electrons have r s < 1, 
while the valence electrons have r 3 > 1, with a tail stretching toward r, — > oc. To see the 
distribution of densities better, we define the density ofr s 's: 

9i(r s ) = jd\ n(r) *(r. - r.(r)). (15.1) 
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r 

Figure 15.2: r 3 (r) in Ar atom. 



This has the simple interpretation: g{r s )dr s is the number of electrons in the system with 
Seitz radius between r s and r s + dr s , and so satisfies: 

/■CO 

f o dr s g(r s ) = N. (15.2) 
This is plotted in Fig. 15.3. The different peaks represent the s electrons in each shell. If 
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Figure 15.3: gi(r s ) in Ar atom. 

we integrate this function forward from r s = 0, we find that it reaches 2 at r s = 0.15, 4 
at 0.41, 10 at 0.7, and 12 at 0.88; these numbers represent (roughly) the maximum r s in a 
given shell. The Is core electrons live with r s < 1.5, and produce the first peak in g\\ the 
2s produce the peak at r s = 0.2, and the 2p produce no peak, but stretch the 2s peak upto 
0.7. The peak at about 0.82 is the 3s electrons, and the long tail includes the 3p electrons. 
Why is this analysis important? We usually write 

E^ A = J d 3 r n(v) C f (r s (r)), (15.3) 
where e™ lf is the exchange-correlation energy per particle in the uniform gas. But armed with 
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our analysis, we may rewrite this as 

Sx L c DA = / °°^5(r s )C f (r s ). (15.4) 

Thus the contribution to for a given r s value is given by gi(r s ) times the weighting 

factor e™ lf (r s ). Clearly, from Fig. 15.3, we see that, to get a good value for the exchange- 
correlation energy of the Ar atom, LDA must do well for r s < 2, but its performance for 
larger r s values is irrelevant. 

Typical r s values are small for core electrons (at the origin, a hydrogenic atom has r s = 
0.72/Z), but valence electrons have r s between about 1 and 6. These produce the dominant 
contribution to chemical processes, such as atomization of molecules, but core relaxations 
with r s <g 1 can also contribute. The valence electrons in simple metal solids have r s between 
2 and 6. 

15.3 Energy density 



15.4 Questions 
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Chapter 16 

Exchange-correlation hole 



1 In this chapter, we explore in depth the reliability of LSD. Although not accurate enough for 
thermochemistry, LSD has proven remarkably systematic in the errors that it makes. Improved 
functionals should incorporate this reliability. We will see that in fact, understanding the 
inherent reliability of LDA is the first step toward useful generalized gradient approximations. 

16.1 Density matrices and holes 

We now dig deeper, to understand better why LSD works so reliably, even for highly inho- 
mogeneous systems. We must first return to the many-body wavefunction. We define the 
following (first-order) density matrix: 

j(x, x') = N j dx 2 ... J dx N x 2 , ■ ■ ■ , x N )ty(x', x 2 , . . . , x N ) (16.1) 

The diagonal elements of this density matrix are just the spin-densities: 

7(2, x) = n(x) (16.2) 

and the exact kinetic energy can be extracted from the spin-summed (a.k.a reduced) density 
matrix: 

7 (r,r') = £7(r<7,rV) (16 . 3) 

and 

T = -^jd 3 rV 2 1 (T,r')\ r=r , (16.4) 

Note that, although 7(2, x') is a one-body property, only its diagonal elements are deter- 
mined by the density. For example, we know that 7 S , the density matrix in the Kohn-Sham 
system, differs from the true density matrix, since the true kinetic energy differs from the 
non-interacting kinetic energy: 

To = -\j dVV 2 ( 7 (r, r') - 7s (r, r')) | r=r < (16.5) 

1 ©2000 by Kieron Burke. All rights reserved. 
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For a single Slater determinant, we have the following simple result: 

ls (x,x') =<WEC«<Mr') (16.6) 

i=l 

where (j>i a (r) is the i-th Kohn-Sham orbital of spin a. This is diagonal in spin. 

Exercise 70 For same-spin electrons in a large box (from to L, where L — > 00) in one- 
dimension, show that 

, m _ kF ( sm(k F (x - x']) _ sm(k F {x + x')) 
' % \ k F (x — x') kp(x + x') 

Prove this density matrix has the right density, and extract t s (x) from it. 

We will, however, focus more on the potential energy, since that appears directly in the 
usual adiabatic connection formula. To this end, we define the pair density: 

P(x,x') = N(N-l) Jdx 3 ...J dx N \^(x,x',x 3 ,...,x N )f (16.8) 

which is also known as the (diagonal) second-order density matrix. This quantity has an 
important physical interpretation: P(ra,r'a')d 3 rd 3 r' is the probability of finding an electron 
of spin (7 in d 3 r around r, and a second electron of spin a' in dV around r'. Thus it contains 
information on the correlations among the electrons. To see this, we may write: 

P(x,x')d 3 rd 3 r' = {n(x)d 3 r} {n 2 (x,x')d 3 r'} (16.9) 

where n 2 (x,x')d 3 r' is the conditional probability of finding the second electron of spin a' 
in d 3 r' around r', given that an electron of spin a has already been found in d 3 r around r. 
Now, if the finding of the second electron were completely independent of the first event, 
then n 2 (x, x') = n(x'). This can never be true for any wavefunction, since, by definition: 

J dx! P(x, x') = (N - 1) n(x) (16.10) 

implying: 

j dx' n 2 (x,x') = N- 1, (16.11) 

i.e., the fact that we have found one electron already, implies that the remaining conditional 
probability density integrates up to one less number of electrons. 

Exercise 71 Show that 

j dx'P(x, x') = (N - 1) n(x) (16.12) 

and state this result in words. What is the pair density of a one-electron system? Deduce 
what you can about the parallel and antiparallel pair density of a spin-unpolarized two electron 
system, like the He atom. 
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Why is this an important quantity? Just as the kinetic energy can be extracted from 
the first-order density matrix, the potential energy of an interacting electronic system is a 
two-body operator, and is known once the pair density is known: 

V ee = ^JdxJ dx'P{x,x')/\r-r'\ (16.13) 

Although it can be very interesting to separate out the parallel and antiparallel spin contribu- 
tions to the pair density, we won't do that here, since the Coulomb repulsion does not. The 
reduced pair density is the spin-summed pair density: 

P(r,r') = £P(r<7,rV), (16 . 14) 
and is often called simply the pair density. Then 

Vm = lj*rJW^l. (16.15) 

As we have seen, a large chunk of V ee is simply the Hartree electrostatic energy U, and 
since this is an explicit density functional, does not need to be approximated in a Kohn-Sham 
calculation. Thus it is convenient to subtract this off. We define the (potential) exchange- 
correlation hole density around an electron at r of spin a: 

P(x, x') = n{x) {n{x) + n xc (x, x')) . (16.16) 

The hole is usually (but not always) negative, and integrates to exactly — 1: 

Jdx' n xc (x;x') = -1 (16.17) 

Since the pair density can never be negative, 

n xc {x,x') > -n(x'). (16.18) 

This is not a very strong restriction, especially for many electrons. Again, while the spin- 
decomposition of the hole is interesting, our energies depend only on the spin-sum, which is 
less trivial for the hole: 

n xc (r,r') = ^>(z)n xc (z, a/)) /n(r). (16.19) 

The pair density is a symmetrical function of r and r': 

P(r',r) = P(r,r'), (16.20) 

but the hole is not. It is often useful to define u = r' — r as the distance away from the 
electron, and consider the hole as function of r and u. These ideas are illustrated in Fig. 
(16.1). It is best to think of the hole as a function of u = x' — x, i.e., distance from the 
electron point. The hole is often (but not always) deepest at the electron point, and decays 
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Figure 16.1: Cartoon of hole in one-dimensional exponential density. Plotted are the actual density n(x') (long dashes), the 
conditional density 712(3; — 2, a:') (solid line), and their difference, the hole density. n xc (x — 2.x') (dashed line). 

rapidly with distance, as the electrons avoid the electron at x. The product of densities in 
P(r,r') gives rise to U, the Hartree energy, while 

Uxc = I d\ n(r) I d 3 u WXC ^' U) . (16.21) 

Thus we may say, in a wavefunction interpretation of density functional theory, the exchange- 
correlation energy is simply the Coulomb interaction between the charge density and its 
surrounding exchange-correlation hole. 

The exchange hole is a special case, and arises from the Kohn-Sham wavefunction. For a 
single Slater determinant, one can show 

P,(x, x') = £ £ C(r)^V(r') (<Mr)<Mr') - <MO<Mr)) (16.22) 

2=1 3=1 

The first (direct) term clearly yields simply the product of the spin-densities, n(x)n(x'), while 
the second (exchange) term can be expressed in terms of the density matrix: 

P x {x, x') = n(x)n(x') - | 7s (x, x')\ 2 (16.23) 

yielding 

n x (x,x') = -\~/ s {x,x')\ 2 /n(x). (16.24) 

This shows that the exchange hole is diagonal in spin (i.e., only like-spins exchange) and is 
everywhere negative. Since exchange arises from a wavefunction, it satisfies the sum-rule, so 
that 

Jd 3 u n x (r,u) = -1 (16.25) 

The exchange hole gives rise to the exchange energy: 

1 i i to. n(r)n x (r, r') 
E x = -Jdh Jd\' [ >*\; (16.26) 



16.2. HOOKE'S ATOM 



125 



Exercise 72 Spin-scaling the hole 
Deduce the spin-scaling relation for the exchange hole. 

Exercise 73 Holes for one-electron systems 

For the H atom (in 3d), plot P(r, r'), n 2 (r,r') and n x (r, r') for several values of r as a 
function ofr', keeping r' parallel to r. 

The correlation hole is everything not in the exchange hole. Since the exchange hole 
satisfies the sum-rule, the correlation hole must integrate to zero: 

Jd 3 u n c (r,u) = 0. (16.27) 

This means the correlation hole has both positive and negative parts, and occasionally the 
sum of exchange and correlation can be positive. It also has a universal cusp at u = 0, due 
to the singularity in the electron-electron repulsion there. 

An alternative way of representing the same information that is often useful is the pair- 
correlation function: 

g{x,x') = P{x,x')/(n(x)n{x')) (16.28) 

This pair-correlation function contains the same information as the hole, but can make some 
results easier to state. For example, at small separations, the Coulomb interaction between 
the electons dominates, leading to the electron-electron cusp in the wavefunction: 

dg s ^(r, u)/du\ u=Q = g^(r, u = 0), (16.29) 

where the superscript denotes a spherical average in u. Similarly, at large separations, g — > 1 
as u — > oc in extended systems, due to the screening effect, but not so for finite systems. 
Furthermore, the approach to unity differs between metals and insulators. 

16.2 Hooke's atom 

A useful and interesting alternative external potential to the Coulomb attraction of the nucleus 
for electrons is the harmonic potential. Hooke's atom consists of two electrons in a harmonic 
well of force constant k, with a Coulomb repulsion. It is useful because it is exactly solvable, 
so that many ideas can be tested, and also (more importantly), many general concepts can 
be illustrated. 

Exercise 74 Hooke's atom: Approximate HF 

Repeat the approximate HF calculation of chapter ?? on Hooke's atom, using appropriate 
orbitals. Make a definite statement about the ground-state and HF energies. 

Exercise 75 Hooke's atom: Separation of variables 

Show that, for two electrons in external potential v ext (r) = kr 2 /2, and with Coulomb inter- 
action, the ground-state wavefunction may be written as 

*(ri,r 2 ) = <b(R)4>(u), (16.30) 
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where R = (ri + r 2 )/2, u = r 2 — r 1; < &(i?) is the ground-state orbital of a 3-d harmonic 
oscillator of mass 2, and <fi(u) satisfies 

l-V 2 u + h f + lU(u) = e u cj>(u) (16.31) 

Exercise 76 Hooke's atom exactly 

Show that the function 4>(u) = C(l + u/2) exp(— u 2 /4) satisfies the Hooke's atom equation 
(16.31) with k = 1/4, and find the exact ground-state energy. How big was your error in 
your estimated energy above? Make a rigorous statement about the correlation energy. 

Exercise 77 Hooke's atom: Exchange and correlation holes 

Use the exact wavefunction for Hooke's atom at k = 1/4 to plot both n x (r, u) and n c (r, u) 
for r = and 1, with u parallel to r. The exact density is: 

n(r) = v / 2C 2 exp(-r 2 /2)/r x 

{7r + r 3 + 8r exp(-r 2 /2)/v / 2^ : + 4(1 + r 2 )erf(r/v / 2)} (16.32) 

where the error function is 

erf(x) = -J= J* dy exp(-y 2 ) (16.33) 

and C- 1 = 2^571 + 8^ 

Exercise 78 Hooke's atom: Electron-electron cusp 

Repeat the above exercise for the pair-correlation function, both exchange and exchange- 
correlation. Where is the electron-electron cusp? What happens as u — > oc? 

16.3 Transferability of holes 

The pair density looks very different from one system to the next. But let us consider two 
totally different systems: the 1-d H-atom and the (same-spin) 1-d uniform electron gas. For 
any one-electron system, the pair density vanishes, so 

n x (r,r') = -n(r') (N = 1). (16.34) 

For our Id H atom, this is just n(x + u), where n(x) = exp(— 2\x\), and is given by the solid 
curve in Fig. 16.2, with x = 0. On the other hand, we can deduce the exchange hole for 
the uniform gas from the bulk value of the first-order density matrix. From Eq.(16.7), taking 
x,x' to be large, we see that the density matrix becomes, in the bulk, 

7s unif ( u ) = n sm(k F u)/(k F u) (16.35) 

leading to 

n x (u) = — n sin 2 (kFu)/(kFu) 2 (16.36) 
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Figure 10. 2: Exhange hole at center of one-dimensional H-atom, both exactly and in LDA. 

These depend only on the separation between points, as the density is constant in the system. 
They also are symmetric in u, as there cannot be any preferred direction. This hole is also 
plotted in Fig. 16.2, for a density of n = 1, the density at the origin in the 1-d H-atom. 

Comparison of the uniform gas hole and the H-atom hole is instructive. They are remarkably 
similar, even though their pair densities are utterly different. If the electron-electron repulsion 
is taken as 5{u), U = 1/8 for the Id H-atom, while U diverges for the uniform gas. Why are 
they so similar? Because they are both holes of some quantum mechanical system. So both 
holes are normalized, and integrate to -1. They are also equal at the ontop value, u = 0. 

Exercise 79 Ontop exchange hole Show that the ontop exchange hole is a local-spin 
density functional, and give an expression for it. 

Why is this important? We may think of the local approximation as an approximation to 
the hole, in the following way: 

n^ DA (x, x + u) = n™ if (n(x); u), (16.37) 

i.e., the local approximation to the hole at any point is the hole of a uniform electron gas, 
whose density is the density at the electron point. Fig 16.2 suggests this will be a pretty good 
approximation. Then the energy per electron due to the hole is just 



with 



Since, for the uniform gas, 



/oc 
du v ee (u) n x (x,x + u), 
-oc 

/oc 
dx n(x) € x (x). 
-oc- 

e » nif („) = r du v ee {u) < nif (n;tt), 

J—OO 



(16.38) 



(16.39) 



(16.40) 



we may consider E]f k [n] as arising from a local approximation to the hole at x. Thus, the 
general similarities tell us that LDA should usually be in the ball park. In fact, if v ee (u) = 
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—S(u) as we have used in the past, and since the exact ontop hole is just — n(x), then LDA 
would be exact for X in that Id world. But since the real Coulomb interaction in 3d is not a 
contact interaction, we will imagine a Id world in which the repulsion is exp(— 2|u|). Then, 
because the overall hole shapes are similar, we find e(0) = —0.5 exactly, and = —0.562 in 
LDA. 





o r 




-0.2 - 


s 


-0.4 - 








X 

e 


-0.6 - 




-0.8 - 




-1 - 




Figure 16.3: Exhange hole at x — 0.5 of one-dimensional H-atom, both exactly and in LDA. 

We have shown above how to understand qualitatively why LDA would give results in the 
right ball-park. But can we use this simple picture to understand quantitatively the LDA 
exchange? If we calculate the exchange energy per electron everywhere in the system, we 
get Fig. 16.4. This figure resembles that of the real He atom given earlier. Looked at more 

o 




Figure 16.4: Exhange energy per electron in the one-dimensional H-atom, both exactly and in LDA, for repulsion exp(— 2\u\). 

closely, our naive hopes are dashed. In fact, our x = case is more an anomaly than typical. 
We find, for our exponential repulsion, E x = —0.351 in LDA, and —0.375 exactly, i.e., the 
usual 10% underestimate in magnitude by LDA. But at x = 0, the LDA overestimates the 
contribution, while for large x, it underestimates it. 

At this point, we do well to notice the difference in detail between the two holes. The 
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exact H-atom hole contains a cusp at u = 0, missing from the uniform gas hole. This causes 
it to deviate quickly from the uniform gas hole near u = 0. At large distances, the H-atom 
hole decays exponentially, because the density does, while the uniform gas hole decays slowly, 
being a power law times an oscillation. These oscillations are the same Friedel oscillations 
we saw in the surface problem. For an integral weighted exp(— 2|it|), the rapid decay of the 
exact hole leads to a smaller energy density than LDA. 

Now watch what happens when we move the hole point off the origin. In Fig. 16.3, we 
plot the two holes at x = 0.5. The H-atom hole is said to be static, since it does not change 
position with x. By plotting it as a function of it, there is a simple shift of origin. The hole 
remains centered on the nucleus, which is now at u = —0.5. On the other hand, we will 
see that for most many electron systems, the hole is typically quite dynamic, and follows the 
position where the first electron was found. This is entirely true in the uniform gas, whose 
hole is always symetrically placed around the electron point, as seen in this figure. Note that, 
although it is still true that both holes are normalized to — 1, and that the ontop values still 
agree, the strong difference in shape leads to more different values of e(x) (-0.312 in LDA, 
and -0.368 exactly). In Fig. 16.4, we plot the resulting exchange energy/electron throughout 
the atom. The LDA curve only loosely resembles the exact curve. Near the nucleus it has a 
cusp, and overestimates e x (x). As x gets even larger, n{x) decays exponentially, so that the 
LDA hole becomes very diffuse, since its length scale is determined by l/k F , where k F = irn, 
while the exact hole never changes, but simply moves further away from the electron position, 
large x. When weighted by the electron density, to deduce the contribution to the exchange 
energy, there is a large cancellation of errors between the region near the nucleus and far 
away. (To see the similarities with real systems, compare this figure with that of Fig. 17.7). 




Figure 16.5: Symmetrized exhange hole at x — 0.5 of one-dimensional H-atom, both exactly and in LDA. 



But, while all the details of the hole are clearly not well-approximated in LDA, especially as 
we move around in a finite system, we now show that the important averages over the hole 
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are. To begin with, both positive and negative values contribute equally, so we may write: 

/CO roc 
dxn(x) dun x Ym (x,u), (16.41) 
-co J0 

where 



n s ^ m (x, u) = \ (n x (x, u) + n x (x, —it)) 



(16.42) 



This symmetrizing has no effect on the LDA hole, since it is already symmetric, but Fig. 
?? shows that it improves the agreement with the exact case: both holes are now parabolic 
around u = 0, and the maximum deviation is much smaller. However, the cusp at the nucleus 
is clearly missing from the LDA hole, and still shows up in the exact hole, 
o 




Figure 10. G: System-averaged symmetrized exhange hole of one-dimensional H-atom, both exactly and in LDA. 

Our last and most important step is to point out that, in fact, it is the system-averaged 
symmetrized hole that appears in E x , i.e., 

/CO 
ao dxn(x)n s x Ym (x;u) (16.43) 

where now it always taken to be positive, since 

roo 

E x = j o du v ee (u)(n x (u)) (16.44) 

The system-averaged symmetrized hole is extremely well-approximated by LDA. Note how no 
cusps remain in the exact hole, because of the system-averaging. Note how both rise smoothly 
from the same ontop value. Note how LDA underestimates the magnitude at moderate it, 
which leads to the characteristic underestimate of LDA exchange energies. Finally, multiplying 
by fee = exp(— 2|it|), we find Fig. 16.7. Now the 10% underestimate is clearly seen, with no 
cancellation of errors throughout the curve. 

We can understand this as follows. For small u, a local approximation can be very accurate, 
as the density cannot be very different at x + u from its value at x. But as it increases, the 
density at x + u could differ greatly from that at x, especially in a highly inhomogeneous 
system. So the short-ranged hole is well-approximated, but the long-range is not. In fact, for 
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Figure 16.7: System-averaged symmetrized exhange hole of one-dimensional H-atom, both exactly and in LDA, weighted by 
exp(— 2\u\). 

large u, the oscillating power-law LDA behavior is completely different from the exponentially 
decaying exact behavior: 

Exercise 80 Show that the exact system-averaged symmetrized hole in the 1-d H atom is: 

(n x (u)) = -exp(-2w) (l + 2u), (= f(iru) LSD), (16.45) 

where 

f( z ) = -y 2 .fo dz ' sinV)/V. (16.46) 
But the constraint of the sum-rule and the wonders of system- and spherical- (in the 3d case) 
averaging lead to a very controlled extrapolation at large u. This is the true explanation of 
LDA's success. Notice also that this explanation requires that uniform gas values be used: 
Nothing else implicitly contains the information about the hole. 



16.4 Old faithful 



We are now in a position to understand why LSD is such a reliable approximation. While it 
may not be accurate enough for most quantum chemical purposes, the errors it makes are 
very systematic, and rarely very large. 

We begin from the ansatz that LSD is a model for the exchange-correlation hole, not just 
the energy density. We denote this hole as n™ lf (r s , u), being the hole as a function of 
separation u of a uniform gas of density (47rrf/3) _1 and relative spin-polarization £. Then, 
by applying the technology of the previous section to the uniform gas, we know the potential 
exchange-correlation energy density is given in terms of this hole: 

<f (r„ C) = 2tt duu nlf (r„ C; «). (16.47) 

This then means that, for an inhomogeneous system, 

E4 S >] = / d 3 r <f (r,(r), C(r)) = N / o °° duu (16.48) 
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KfM) = ^Jdh n(r) ^<f(r s (r),C(r); M ) 



N J " ' " v_/ 4tt 
is the system- and spherically-averaged hole within LSD. 



(16.49) 
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Figure 16.8: Universal curve for the system-averaged on-top hole density in spin-unpolarized systems. The solid curve is for 
the uniform gas. The circles indicate values calculated within LSD, while the crosses indicate essentially exact results, and 
the plus signs indicate less accurate CI results. 

Why should this hole look similar to the true hole? Firstly, note that, because the uniform 
gas is an interacting many-electron system, its hole satisfies the same conditions all holes 
satisfy. Thus the uniform gas exchange hole integrates to -1, and its correlation hole integrates 
to zero. Furthermore, its exchange hole can never be negative. Also, the ontop hole (u = 0) 
is very well approximated within LSD. For example, in exchange, 



n x (r,r) 



f(r) + nf(r))/n(r) 



(16.50) 



i.e., the ontop exchange hole is exact in LSD. Furthermore, for any fully spin-polarized or 
highly-correlated system, 

n xc (r, r) = -n(r) (16.51) 

i.e., the ontop hole becomes as deep as possible (this makes P(r, r) = 0), again making 
it exact in LSD. It has been shown to be highly accurate for exchange-correlation for most 
systems, although not exact in general. Then, with an accurate ontop value, the cusp 
condition, which is also satisfied by the uniform gas, implies that the first derivative w.r.t. u 
at u = is also highly accurate. In Fig. 16.8, we plot the system-averaged ontop hole for 
several systems where it is accurately known. For these purposes, it is useful to define the 
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system -averaged density: 

<«(«)> = ^ 

and a system-averaged mean r, value 



47T 



y <fV n(r) n(r + u). 



/ d 3 r n 2 (r) r s (r) 



(16.52) 



(16.53) 



/ c( 3 r n 2 (r) 

Clearly, the system-averaged ontop hole is very accurate in LSD. 

That LSD is most accurate near u = can be easily understood physically. At u grows r' 
gets further and further away from r. But in LSD, our only inputs are the spin-densities at r, 
and so our ability to make an accurate estimate using only this information suffers. Indeed, as 
mentioned above, at large separations the pair-correlation function has qualitatively different 
behavior in different systems, so that LSD is completely incorrect for this quantity. Thus, we 
may expect LSD to be most accurate for small u, as indeed it is. But even at large u, its 
behavior is constained by the hole normalization sum-rule, 
o 




Figure 16.9: System-averaged exchange hole density (in atomic units) in the He atom, in LDA (dashed) and exactly (HF - 
solid) . The area under each curve is the exchange energy. 

Two other points are salient. The exchange-correlation hole in the uniform gas must be 
spherical, by symmetry, whereas the true hole is often highly aspherical. But this is irrelevant, 
since it is only the spherical-average that occurs in E xc . Furthermore, the accuracy of LSD 
can fail in regions of extreme gradient, such as near a nucleus or in the tail of a density. But 
in the former, the phase-space in the system-average is small, while in the latter, the density 
itself is exponentially small. The same argument applies to large separations: even if g is 
badly approximated by LSD, it is — n(r + u)g(r, r + u) that appears in the hole, and this 
vanishes rapidly, both exactly and in LSD. Thus limitations of LSD in extreme situations do 
not contribute strongly to the exchange-correlation energy. 

This then is the explanation of LSD's reliability. The energy-density of the uniform gas 
contains, embedded in it, both the sum-rule and the accurate ontop hole information. For 
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Figure 16.10: System-averaged correlation hole density at full coupling strength (in atomic units) in the He atom, in LSD, 
numerical GGA, and exactly (CI). The area under each curve is the full coupling-strength correlation energy. 

most systems, this implies that their exact system-averaged hole looks very like its LSD 
approximation, with only small differences in details, but never with very large differences. 
In Figs. 16.9-16.10, we plot the system-averaged exchange and potential-correlation holes 
for the He atom, both exactly and in LSD. Each is multiplied by a factor of 2nu, so that 
its area is just the exchange or correlation energy contribution (per electron). We see that 
LSD is exact for small u exchange, and typically underestimates the hole in the all-important 
moderate u region, while finally dying off too slowly at large u. For correlation, we see that 
the LSD hole is always negative in the energetically significant regions, whereas the true hole 
has a significant positive bump near u = 1.5. This is yet another explanation of the large 
overestimate of LSD correlation. In the next chapter, we will show how to use the failed 
gradient expansion to improve the description of the hole, and so construct a generalized 
gradient approximation. 



Part IV 
Beyond LDA 
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Gradients 



In this chapter, we explore the "obvious" correction to the local approximation, namely 
the inclusion of gradient information. We '11 see that its not as easy as it looks, and why 
it took about a quarter of a century to accomplish. 

17.1 Perimeter problem 

To illustrate the interesting complexities of the use of gradients, we return to our problem 
from Chapter 2, in which we needed to know the perimeter of a shape, given its definition 
in terms of polar coordinates, r{6). In all of what follows, we assume we do not know the 
exact answer. 

We already saw in chapter 8 that the local approximation is 

P loc = 06 r{ff). (17.1) 
Let us first consider the family of smooth curves introduced in that chapter: 

r(B) = 1 + e cos(n#) (17.2) 

Here < e < 1, while n is an integer. For these curves, the local approximation yields exactly 
27r for the perimeter. 

In table 17.1, I have compared results for n = 1, as a function of e. The local approximation 
is doing rather well, with a maximum error of —15%. Notice how its always an underestimate. 



n 


e 


local 


% error 


exact 


1 


0.0 


6.2832 


0.00 


6.2832 


1 


0.1 


6.2832 


-0.25 


6.2989 


1 


0.2 


6.2832 


-0.99 


6.3462 


1 


0.4 


6.2832 


-3.88 


6.5371 


1 


0.8 


6.2832 


-14.37 


7.3376 



Table 17.1: Perimeter's of shapes for n — 1, and various e, both exactly and in local approximation. 
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0.2 


6.2832 


-1 


6.3462 
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0.2 


6.2832 


-4 


6.5297 


4 


0.2 


6.2832 


-13 


7.1989 


8 


0.2 


6.2832 


-32 


9.2984 


16 


0.2 


6.2832 


-57 


14.7104 


32 


0.2 


6.2832 


-77 


26.7835 


64 


0.2 


6.2832 


-88 


51.9081 



Table 17.2: Perimeter's of shapes for e — 0.2, and various n, both exactly and in local approximation. 

One could imagine a clever person, not knowing the true formula, being able to prove an 
exact condition like: 

P loc < P (17.3) 

We also notice that the error increases with e. Obviously e = is a circle, where the local 
approximation is exact. For small e, the curve is in some sense close to a circle, and the 
local approximation should be accurate. In Fig. 17.1, we plot the e = 0.8 case, and see that, 

|n=1 , eps=0.8 I 



Figure 17.1: Shape r = 1 + O.8cos(0). 

although the radius varies from 0.2 to 1.8, the curve looks quite circular. 

On the other hand, in Table 17.1, we fix e = 0.2, and increase n. The error grows almost 
quadratically with n. Again, n = is the circular case, and so small n is somehow more 
circular than large n. 

To understand the ever increasing error, we plot the n = 64 curve in Fig. 17.2. The local 
approximation yields the perimeter of a circle with the average radius, r = 1. Obviously, the 
64 wiggles between r = 1.2 and r = 0.8 as one goes round the circle are not accounted for 
in the local approximation, but add greatly to the perimeter. 

Thus we can speak of two distinct ways in which a curve might be 'close' to a circle. The 
first, more familiar way, is when e is small, and so the perturbation on r = 1 is weak. This can 
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Figure 17.2: Shape r = 1 + 0.2cos(646>). 

usually be treated by response theory. Perturbation theory will always be accurate once the 
pertubation is sufficiently week, no matter how rapidly the curve is varying (i.e., no matter 
how large n is). The second, less familiar way is when the perturbation is slowly-varying, but 
can be arbitrarily large. This is the case when n is small, but e need not be. In Fig. 17.1, 
the local approximation works quite well, even though e is enormous. The ratio of maximum 
to minimum r is 9! 

Let us first analyze the weak perturbation case. We can cheat by using our knowledge of the 
exact functional, but I'm sure there are ways to derive the result without it. If r = r + e f(6), 
then 

p = 2 ^+lU 2 MS) 2+0(e4) (l7 - 4) 

For our standard curve, the integral is r?, 2 /T. There is no linear term, because it vanishes by 
periodicity requirements. You will find this formula quite accurate, even at e = 0.2, and 
hence the (near) quadratic growth in the error in the local approximation in Table 17.1. 

Note that, in the local approximation, there is no e 2 term. Thus the local approximation, 
while being exact for the circle, is hopeless for weak perturbations around the circle. 

The other good approximation is called the gradient expansion. We first note that the 
only way to make a dimensionless gradient is by dividing the derivative by the radius. We 
define s = dr/dO/r. Then, for a slowly-varying shape, we can write 

P = J d6r(l + Cs 2 + ...) (17.5) 

where C is yet to be determined. There is no term linear in s, because its integral would 
vanish. Thus the GEA, or gradient expansion approximation, consists of keeping just the first 
two terms, once C has been found. 

Finding C is easy, once we know the linear response. We simply imagine a perturbation 
that is both weak and slow. Then s = dr/d9/r Q , and we see that C must be 1/2. Thus the 
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Table 17.3: Perimeter's of shapes for e = 0.2, and various n, both exactly and in local approximation and in GEA. 

GEA in this case is 

P^P- + i/>i(§) ! ,1,6, 

In Table 17.1, we have added the results of the gradient expansion. We see that for 
n = 1, where the local approximation was quite good, GEA reduces the error by at least 
a factor of 4, and sometimes much more. It also overestimates the perimeter in all cases. 
Even up to n = 8, for e = 0.2, the error is still reduced by a factor of 3. But for larger n, 
meaning larger gradients, the GEA overcorrects, eventually producing larger errors than the 
local approximation! Our conclusion is that, for slowly-varying shapes, the gradient expansion 
greatly improves on the local approximation, but does not work for rapid variations, and can 
even worsen the results. 

To make the point very clear, suppose we live in a world where most of the shapes we care 
about have n between 2 and 8, while e varies between .4 and .8. The results of the local and 
gradient expansion approximations are listed in Table ??. Only for the most slowly-varying 
cases does the GEA really improve over the local approximation. In all cases, it overcorrects, 
usually by more than the original error. For these systems, the GEA is hardly an improvement. 

17.2 Gradient expansion 

Way back when, in the original Kohn-Sham paper, it was feared that LSD might not be too 
good an approximation (it turned out to be one of the most successful ever), and a simple 
suggestion was made to improve upon its accuracy. The idea was that, for any sufficiently 
slowly varying density, an expansion of a functional in gradients should be of ever increasing 
accuracy: 

A GEA [n] = J d 3 r (o(n(r)) + &(n(r))| Vn| 2 + . . .) (17.7) 
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Table 17.4: Perimeter's of shapes for various e and n, both exactly and in local approximation and in GEA. 

Then, if LSD was moderately accurate for inhomogeneous systems, GEA, the gradient ex- 
pansion approximation, should be more accurate. The form of these gradient corrections can 
be determined by scaling, while coefficients can be determined by several techniques. 

For both T s and E x , defined on the Kohn-Sham wavefunction, the appropriate measure of 
the density gradient is given by 

s(r) = |Vn(r)|/(2fc F (i>(r)) (17.8) 

This measures the gradient of the density on the length scale of the density itself, and has 
the important property that s[n 7 ](r) = s[n](7r), i.e., it is scale invariant. It often appears in 
the chemistry literature as x = | Vn|/(n 4 / 3 ), so that x « 6s. The gradient expansion of any 
functional with power-law scaling may be written as: 

A GEA [n] = J d 3 r a(n(r)) (l + Cs 2 (r)) (17.9) 

The coefficient may be determined by a semiclassical expansion of the Kohn-Sham density 
matrix, whose terms are equivalent to a gradient expansion. One finds 

C s = 5/27, C x = 10/81 (17.10) 

(In fact, a naive zero-temperature expansion gives C x = 7/81, but a more sophisticated 
calculation gives the accepted answer above). The gradient correction generally improves 
both the non-interacting kinetic energy and the exchange energy of atoms. The improvement 
in the exchange energy is to reduce the error by about a factor of 2. 

Exercise 81 For the H-atom, calculate the gradient corrections to the kinetic and exchange 
energies. (Don't forget to spin-scale). 
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For correlation, the procedure remains as simple in principle, but more difficult in practice. 
Now, there are non-trivial density- and spin- dependences, and the terms are much harder to 
calculate in perturbation theory. In the high-density limit, the result may be written as 

AE ° EA = 3^/ dh ™( r MCM)i 2 (r)- (17.11) 

Here t = \Vn\/(2k s n), where k s = ^jAkpfii is the Thomas-Fermi screening length, the 
natural wavevector scale for correlations in a uniform system. The spin-polarization factor is 
<KC) = ((1 + C) 2,/3 + (1 — C) 2 ^ 3 ) /2- This correction behaves poorly for atoms, being positive 
and sometimes larger in magnitude than the LSD correlation energy, leading to positive 
correlation energies! 

Exercise 82 GEA correlation energy in H atom 

Calculate the GEA correction to the H atom energy, and show how AE GEA scales in the 
high-density limit. 

17.3 Gradient analysis 

Having seen in intimate detail which r s values are important to which electrons in Fig. 15.2, 
we next consider reduced gradients. Fig. 17.3 is a picture of s(r) for the Ar atom, showing 
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Figure 17.3: s(r) in Ar atom. 

how s changes from shell to shell. Note first that, at the nucleus, s has a finite moderate 
value. Even though the density is large in this region, the reduced gradient is reasonable. 

Exercise 83 Reduced gradient at the origin 

Show that at the origin of a hydrogen atom, s(0) = 0.38, and argue that this value wont 
change much for any atom. 

For any given shell, e.g. the core electrons, s grows exponentially. When another shell begins 
to appear, there exists a turnover region, in which the gradient drops rapidly, before being 
dominated once again by the decay of the new shell. We call these intershell regions. 
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We define a density of reduced gradients, analogous to our density of r s 's: 

g 3 (s) = fd\n{v) 5{s-s{v)), (17.12) 
also normalized to N , and plotted for the Ar atom in Fig. 17.4. Note there is no contribution 
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Figure 17.4: g^(s) in Ar atom. 



to atoms (or molecules) from regions of s = 0. Note also that almost all the density has 
reduced gradients less than about 1.2, with a long tail stretching out of larger s. 

Exercise 84 Gradient analysis for H atom 

Calculate r a (r) and s(r) for the H atom, and make the corresponding gi(r s ) and §3(3) plots. 
Why does the exponential blow up at large r not make the LDA or GEA energies diverge? 

We now divide space around atoms into spheres. We will (somewhat arbitrarily) denote 
the valence region as all r values beyond the last minimum in s(r), containing 6.91 electrons 
in Ar. We will denote the region of all r smaller than the position of the last peak in s as 
being the core region (9.34 electrons). The last region in which s decreases with r is denoted 
the core-valence region, in which the transition occurs (1.75 electrons). We further denote 
the tail region as being that part of the valence region where s is greater than its maximum 
in the core (1.54 electrons). 

To see how big s can be in a typical system that undergoes chemical reactions, consider 
Fig. 17.5. On the left is plotted values of r„ (top) and of s (bottom) as a function of r in an N 
atom. Then on the right are plotted the distribution of values (turned sideways). We see that 
the core electrons have r s < 0.5, while the valence electrons have 0.5 < r s < 2. The s-values 
are more complex, since they are similar in both the core and valence regions. We see that 
almost all the density has 0.2 < s < 1.5. Because s is not monotonic, many regions in space 
can contribute to 53(5) for a given s value. To better distinguish them, one can perform a 
pseudopotential calculation, that replaces the core electrons by an effective potential designed 
only to reproduce the correct valence density. This is marked in the two right-hand panels 
by a long-dashed line. In the 53 (r s ) figure, we see indeed that the pseudopotential is missing 
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Figure 17.5: Distribution of densities and gradients in the N atom. 

precisely the inner core peak, while in the g(s) figure, we see the contribution of the valence 
electrons alone to the peak from the combined core-valence region. 

We may use the s curve to make more precise (but arbitrary) definitions of the various 
shells. Up to the first maximum in s we call the core. The small region where s drops rapidly 
is called the core-valence intershell region. The region in which s grows again, but remains 
below its previous maximum, we call the valence region. The remainer is denoted the tail. 
We make these regions by vertical dashed lines, which have been extended to the upper r s 
curve. Because that curve is monotonic, we can draw horizontal lines to meet the vertical 
lines along the curve, and thus divide the g(r s ) curve into corresponding regions of space. 

But we are more interested in atomization energies than atomic energies. In the case of a 
molecule, there is no simple monotonic function, but one can easily extract the densities of 
variables. Having done this for the molecule, we can make a plot of the differences between 
(twice) the atomic curves and the molecular curve. Now the curves are both positive and 
negative. We see immediately that the core does not change on atomization of the molecule, 
and so is irrelevant to the atomization energy. Interestingly, it is also clear that it is not 
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Figure 17.6: Differences in distribution of densities and gradients upon atomization of the N2 molecule. 

only the valence electrons (r s > 0.8) that contribute, but also the core-valence region, which 
include r s values down to 0.5. Finally, note that the r s values contributing to this energy also 
extend to higher r s regions. 

To summarize what we have learned from the gradient analysis: 

• Note again that nowhere in either the N or Ar atoms (or any other atom) has s << 1, 
the naive requirement for the validity of LSD. Real systems are not slowly-varying. Even 
for molecules and solids, only a very small contribution to the energy comes from regions 
with s < 0.1. 

• A generalized gradient approximation, that depends also on values of s, need only do 
well for s < 1.5 to reproduce the energy of the N atom. To get most chemical reactions 
right, it need do well only for s < 3. 

We end by noting that while this analysis can tell us which values of r s and s are relevant 
to real systems, it does not tell us how to improve on LSD. To demonstrate this, in Fig. 
17.7, we plot the exchange energy per electron in the He atom. One clearly sees an apparent 
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Figure 17.7: Exchange energy per electron in the He atom, both exactly and in LDA. The nucleus is at r — 0. 

cancellation of errors throughout the system, i.e., LDA is not describing this quantity well at 
each point in the system, but it does a good job for its integral. As we will show below, a 
pointwise analysis of energy densities turns out to be, well, somewhat pointless... 



17.4 Questions 

1. If GEA overcorrects LDA, what can you say about how close the system is to uniform? 

2. Do you expect LDA to yield the correct linear response of a uniform electron gas? 

3. Sketch the gradient and density analysis for the ionization energy for Li. 



Chapter 18 

Generalized gradient approximation 



1 In this chapter, we first discuss in more detail why the gradient expansion fails for finite 
systems, in terms of the exchange-correlation hole. Fixing this leads directly to a generalized 
gradient approximation (GGA) that is numerically defined, and corrects many of the limita- 
tions of GEA for finite systems. We then discuss how to represent some of the many choices 
of forms for GGA, using the enhancement factor. We can understand the enhancement factor 
in terms of the holes, and then understand consequences of gradient corrections for chemical 
and physical systems. Finally we mention other popular GGAs not constructed in this fashion. 

18.1 Fixing holes 

We begin with exchange. We have seen that LSD works by producing a good approximation 
to the system- and spherically-averaged hole near u = 0, and then being a controlled extrap- 
olation into the large u region. The control comes from the fact that we are modelling the 
hole by that of another physical system, the uniform gas. Thus the LSD exchange hole is 
normalized, and everywhere negative. When we add gradient corrections, the model for the 
hole is no longer that of any physical system. In particular, its behavior at large separations 
goes bad. In fact, the gradient correction to the exchange hole can only be normalized with 
the help of a convergence factor. 

To see this, consider Fig. 18.1, which shows the spherically averaged exchange hole for a 
given reduced gradient of s = 1. Everything is plotted in terms of dimensionless separation, 
z = 2kpu. The region shown is the distance to the first maximum in the LDA hole, where it 
just touches zero for the first time. The GEA correction correctly deepens the hole at small 
u, but then contains large oscillations at large u. These oscillations cause the GEA hole to 
unphysically become positive after z = 6. These oscillations are so strong that they require 
some damping to make them converge. Without a damping factor, the gradient correction 
to the exchange energy from this hole is undefined. 

The real-space cutoff construction of the generalized gradient approximation (GGA), is 
to include only those contributions from the GEA exchange hole that are negative, and to 

1 ©2000 by Kieron Burke. All rights reserved. 
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Figure 18.1: Spherically-averaged exchange hole density tis PLdv '(u) for s — 1. 

truncate the resulting hole at the first value of u which satisfies the sum-rule. In the figure, for 
small z, the GGA hole becomes more negative than GEA, because in some directions, the GEA 
hole has become positive, and these regions are simply sliced out of GGA, leading to a more 
negative spherically-averaged hole. Then, at about 2k F u = 6.4, the hole is truncated, because 
here the normalization integral equals -1, and its energy density contribution calculated. 
We note the following important points: 

• Simply throwing away positive contributions to the GEA exchange hole looks very ugly 
in real space. But remember we are only trying to construct a model for the spherical- 
and system-average. Thus although that has happened in some directions for small z, 
the spherically-averaged result remains smooth. 

• By throwing out postive GEA contributions, a significantly larger energy density is found 
for moderate values of s. Even for very small s, one finds corrections. Even when GEA 
corrections to the LSD hole are small, there are always points at which the LSD hole 
vanishes. Near these points, GEA can make a postive correction, and so need fixing. 
This means that even as s — » 0, the GGA energy differs from GEA. 

• For very large s (e.g. greater than 3), the wild GEA hole produces huge corrections 
to LSD, the normalization cutoff becomes small, and limits the growth of the energy 
density. However, this construction clearly cannot be trusted in this limit, since even a 
small distance away from u = may have a completely different density. 

The results of the procedure are shown in Fig. 18.2. We see that the LDA underestimate 
is largely cured by the procedure, resulting in a system-averaged hole that matches the exact 
one better almost everywhere. 
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Figure 18.2: System-averaged exchange hole density (in atomic units) in the He atom. The exact curve (solid) is Hartree- 
Fock, the long dashes denote LDA, and the short dashes are real-space cutoff GGA. The area under each curve is the exchange 
energy. 

There is a simple way to picture the results. We have seen how the spin-scaling re- 
lation means that one need only devise a total density functional, and its corresponding 
spin-dependence follows. Furthermore, since it scales linearly, the only way in which it can 
depend on the first-order gradient is through s, the reduced dimensionless gradient. Thus we 
may write any GGA for exchange, that satisfies the linear scaling relation, as 

- unif (n(r))F x (s(r)) (18.1) 



E? GA \n] 
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The factor F x is the exchange enhancement factor, and tells you how much exchange is 
enhanced over its LDA value. It is the analog of F™ f (r s ) from before. In Fig. 18.3, we plot 

1.8 




0.5 1 1.5 2 2.5 3 



s = \Vn\/(2k F n) 

Figure 18.3: Exchange enhancement factors for different GGA's 

various enhancement factors. First note that F x = 1 corresponds to LDA. Then the dashed 
line is GEA, which has an indefinite parabolic rise. The solid line is the result of the real- 
space cutoff procedure. It clearly produces a curve whose enhancement is about double that 
of GEA for small and moderate s, because of the elimination of positive contributions to the 
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hole. But the hole-GGA curve also grows less rapidly at large s, because of the normalization 
cutoff. 

Also included in the figure are two modern GGAs, (Becke 88 and Perdew-Wang 91) that 
we discuss more later. They both agree with the general shape of the hole-GGA for moderate 
values of s, making the real-space cutoff procedure a justification for any of them. The most 
significant deviation in this region is that the hole-GGA has an upward bump near s = 1. 
This is due to the negativity and normalization cutoffs meeting, as is about to happen in Fig. 
?? for slightly larger s, and is an artifact of the crude construction. Recalling that GGA is 
still an approximation, the differences between different GGA's are probably of the order of 
the intrinsic error in this approximate form. 
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Figure 18.4: Spherically-averaged correlation hole density n c for r s = 2 and £ — 0. GEA holes are shown for four values of 
the reduced density gradient, t — |Vn|/(2fc s n). The vertical lines indicate where the numerical GGA cuts off the GEA hole 
to make / Q c dv 47T v 2 n c (v) — 0. 

The case for correlation is very similar. Here there is no negativity constraint, and the 
real-space construction of the GGA correlation hole simply corrects the lack of normalization 
in the GEA correlation hole. Fig. 18.4 shows what happens. The line marked is the LDA 
correlation hole. Note its diffuseness. Here the natural length scale is not kpu, but rather k s u, 
because we are talking about correlation. Similarly reduced gradients are measured relative 
to this length scale. The other curves are GEA holes for increasing reduced gradients. Note 
that the GEA hole is everywhere positive, leading to consistently positive corrections to LDA 
energies. Note also how large it becomes at say t = 1.5, producing a stongly positive energy 
contribution. For small reduced gradients, the gradient correction to the hole is slight, and 
has little effect on the LSD hole, until at large distances, it must be cutoff. As the gradient 
grows, the impact on the LSD hole becomes much greater, and the cutoff shrinks toward 
zero. Since the energy density contribution is weighted toward u = relative to this picture, 
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and all curves initially drop, the GGA correlation hole always produces a negative correlation 
energy density. 
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Figure 18.5: System-averaged correlation hole density at full coupling strength (in atomic units) in the He atom. The exact 
curve (solid) is from a CI calculation, the long dashed is LDA, and the short-dashes are real-space cutoff GGA. The area 
under each curve is the potential contribution to correlation energy. 

In Fig. 18.5, we see the results for the system-averaged hole, which are even more dramatic 
than the exchange case. At small separations, the GEA correclty makes the LDA hole less 
deep, as expected. But now, instead of producing a wild (i.e., unnormalized) postive peak, 
the postive contribution is much less pronounced, and is controlled by the normalization. 
The upward bump and rapid turnoff beyond u = 2 are artifacts of the crude real-space cutoff 
procedure, and remind us of how little we know of the details of the corrections to LDA. But 
the overall effect on the energy is very satisfying. The LDA hole spreads out enormously, 
giving rise to a correlation energy that is at least a factor of 2 too deep, while the real-space 
cutoff hole has produced something similar to the true hole. 

18.2 Visualizing and understanding gradient corrections 

To understand gradient corrections for exchange-correlation, we begin with the spin-unpolarized 
case, where we write 

E^{n] = J cf'r eT il (n(v)) F xc (r s (r), s(r)). (18.2) 
This enhancement factor contains all our others as special cases: 

F xc (r s , 8 = 0) = F»f(r s ), F xc (r s = 0,a) = F x (s). (18.3) 

We plot it for the case of the PBE functional in Fig 18.6. We will discuss the various kinds of 
GGA functionals that have been developed and are in use later, but for now we simply note 
that this GGA was designed to recover the GGA for the real-space correction of the GEA hole 
for moderate s, and includes also several other energetically significant constraints. 
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Figure 18.0: Exchange-correlation enhancement factors for the PBE GGA. 

These curves contain all the physics (and ultimately) chemistry behind GGA's. We make 
: following observations: 

» Along the j/-axis, we have F™ lf (r s ), the uniform gas enhancement factor. 

» The effect of gradients is to enhance exchange. We have already seen this in action in 
the previous sections. Essentially, in the positive direction of the gradient, the density 
is increased, while dropping on the opposite side. This allows the center of the hole to 
move in that direction, becoming deeper due to the higher density, and yielding a higher 
exchange energy density. 

• The effect of gradients is turn off correlation relative to exchange. This is because, in 
regions of high gradient the exchange effect keeps electrons apart, so that their correlation 
energy becomes relatively smaller. In our real-space analysis of the hole, we see that large 
gradients ultimately eliminate correlation. 

» A simple picture of the correlation effects is given by the observation that for t up to 
about 0.5, the GGA make little correction to LDA, but then promptly cuts off correlation 
as t grows, as in Fig. ??. But t is the reduced gradient for correlation, defined in terms 
of the screening wavevector, not the Fermi wavevector. For an unpolarized system: 
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(18.4) 



Thus, for r s about 2, t = s, and we see correlation turning off at t about 0.5. For 
r s = 10, t « 0.4s, so that correlation turns off about s = 1.2, while for r s = 1/2, 
t « 1.7s, so that correlation turns off very quickly, at s about 0.3. 

» An interesting consequence of the point above is that pure exchange is the least local 
curve, and that as correlation is turned on, the functional becomes more local, i.e., closer 
to the LDA value, both in the sense that gradient corrections become significant at larger 
s values, and also that even when they do, their magnitude is less. 
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We can repeat this analysis to understand the effects of spin polarization. Must fill this in. 

Exercise 85 Deduce the maximum value for F xc allowed by the Lieb-Oxford bound. 

Exercise 86 The lines of different r s values never cross in Fig. Y. What does this imply 
about the functional? 

18.3 Effects of gradient corrections 

In this section, we survey some of the many properties that have been calculated in DFT 
electronic structure calculations, and try to understand the errors made LSD, and why GGA 
corrects them the way it does. 

• Both the large underestimate in the magnitude of total exchange energies and the over- 
estimate in the magnitude of correlation energies can be seen immediately from Fig. 
18.6. Ignoring gradients ignores the enhancement of exchange and the truncation of 
correlation. 

• In studying ionization energies of atoms, one finds LSD overbinds s electrons relative to 
p, p relative to d, and so forth. Since core electrons have higher density than outer ones, 
the local approximation is less accurate for them. 

• When a molecule is formed from atoms, the density becomes more homogeneous than 
before. The region of the chemical bond is flatter than the corresponding regions of 
exponential decay in atoms. Thus LDA makes less of an underestimate of the energy of 
a molecule than it does of the constituent atoms. Write 

E xc (system) = E^ A {system) + AE x f A (system) (18.5) 

Then AE xc (atoms) < AE xc {molecules) , remebering that both are negative. The 
atomization energy of a molecule is 

£atmiz = E xc (molecule) - E xc (atoms) (18.6) 

so that 

E xc (atom) = E x ^ A (atom) + AE X ® A (molecule) - E x ^ A (atom) (18.7) 
Thus molecules are overbound in LDA, typically by as much as 30 kcal/mol. 

• Transition state barriers can be understood in much the same way. A typical transition 
state in quantum chemistry is one of higher symmetry and coordination than the reac- 
tants. Thus comparing the transition state to the reactants, in just the way done above 
for a bond and its constituent atoms, we find that the transition-state barrier, defined as 

^trans = _ e xc (re act ants) + E xc (transitionstate) (18.8) 
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will be too low in an LDA calculation, and be raised by a GGA calculation. Notice though 
that this depends on the coordination argument. If the transition state is less coordinated, 
the result is the reverse: LDA overestimates the barrier, and GGA will weaken it. This is 
often the case for, eg. surface diffusion, where the transition state (e.g. bridge-bonded) 
has lower symmetry than the initial state (e.g. atop bonded). 

• By the same reasoning as for atomization energies, consider the process of stretching a 
bond. As a bond is stretched, LDA's overbinding tendancy will reduce. Thus equilibrium 
bond lengths are usually too small in LDA, and GGA stretches them. The only exception 
to this is the case of bonds including H atoms. One can show that typically GGA favors 
a process in which 

f^i^*». <*»> 

where (s), etc., are precisely chosen averages, and where P is typically close to 1, and 
Q w 0. Usually fractional changes in the density are very small, due to the presence of 
core electrons, and the right hand side is thought of as zero. This leads to the generic 
claim that GGA prefers inhomogeneity, i.e., LDA overfavors homogeneity. For example, in 
our cases above, when a molecule is atomized, its mean gradient increases. The inequality 
means that GGA will like this more than LDA, and so have a smaller atomization energy. 
Neglecting the right-hand-side works for all atomization energies, in which the changes 
are large. But for a bond-stretch, the fractional gradient changes are much smaller. In 
the case of H atoms, the density is sufficiently low as to make its fractional change larger 
than that of the gradient, when stretching a bond with an H atom. Thus in that case, 
the gradient corrections have the reverse effect, and those bonds are shortened. 



18.4 Satisfaction of exact conditions 

We can now look back on the previous chapters, and check that the real-space cutoff GGA 
satisfies conditions that LDA satisfies, and maybe a few more. 

• Size-consistency Obviously, GGA remains size consistent. 

• The Lieb-Oxford bound. 

The LO bound gives us a maximum value for the enhancement factor of 1.804. The 
real-space cutoff will eventually violate this bound, but only at large s, where we do not 
trust it anyhow. 

• Coordinate scaling 

For exchange, clear cos only s dependence. 
Show F xc lines never cross. 

• Virial theorem Trivially true. 
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• High-density correlation 

The numerical GGA satisfies this condition. As we show below, many chemical systems 
are close to this, and, by truncating the long-ranged Coulomb correlation hole, numerical 
GGA yields a finite value. 

• Self-interaction error 

Simple density functionals like LDA and GEA naturally have a self-interaction error, 
because they cannot tell when there is only one electron in the system. While GGA 
might improve numerically the value of E x or E c for that one electron, it will still have 
a residual self-interaction error. 

• Symmetry dilemma 

Just as above, nothing in GGAs construction helps here. 

• Potentials 

From the naive sense... 

• Koopman's theorem No real improvement here. 

18.5 A brief history of GGA's 

We have seen how the real-space cutoff construction imposed on the gradient expansion 
produces a numerically-defined GGA. This is not unique, in that, e.g. a smoother cutoff will 
lead to slightly different curves for moderate gradients, and perhaps much different curves 
for large gradients. 

18.6 Questions about generalized gradient approximations 



Chapter 19 
Hybrids 



19.1 Static, strong, and strict correlation 

19.2 Mixing exact exchange with GGA 

19.3 Questions about adiabatic connection formula and hybrids 

1. Does the Hellmann-Feynman theorem apply to the problem of the particle in a potential 

V(x) = exp(— |x|)? 

2. What happens to a system as A — > oc? 

3. What is the formula for i>g Xt [n](r) in terms of a scaled density? 
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Chapter 20 

Orbital functionals 

20.1 Self-interaction corrections 

20.2 Optimized potential method 

20.3 Gorling-Levy perturbation theory 

20.4 Meta GGA's 

20.5 Jacob's ladder 
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Chapter 21 
Time-dependence 



21.1 Schrodinger equation 

We learn in kindergarten that the evolution of the wavefunction is governed by the time- 
dependent Schrodinger equation: 

ffl(t) = H(t) V(t), *(0) given (21.1) 

where H(t) = T + V(t), and a dot implies differentiation with respect to time. We will 
generally be concerned with a time-dependence of the form 

H(t) = H + V(t) (21.2) 

where the unperturbed Hamiltonian, H has eigenstates 

H„\j) = Ej\j) (21.3) 

and V(t) is zero for t < 0. We may then expand the time-dependent wavefunction in terms 
of these unperturbed eigenstates: 

*(*) = £Cj(t) exp(-i^t) \j) (21.4) 

3 

where the phase factor for undisturbed evolution has been extracted. Insertion into Eq. (21.1) 
yields 

* Cj(t) =T,C k exp(iuj jk t) V jk (t), (21.5) 
k 

a coupled set of first-order differential equations, with Vjk = (j\V\k), ujjk = E ] — Ek and 
with initial values 

Cj(0) = (*(0)|j>. (21.6) 

We usually begin in the ground-state, so that C,(0) = Sj . 

To have a feeling for what time-dependent quantum mechanics really means, we immedi- 
ately show some examples for a simple case. We take our particle in a box and subject it to 
a particular driving potential, 

V(xt) = x £{t) (21.7) 
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This is the potential due to an external time-dependent electric field in the dipole approxi- 
mation. 

Exercise 87 Matrix elements of time-dependent potentials 

Calculate the matrix elements Vjk(t) for (a) a harmonic oscillator and (b) a particle in a box, 
when the perturbation is an electric field. 

We choose 

£{t) = A sin(wi) (1 - exp(-wt)) (21.8) 

This form allows a slow turn-on of a periodic electric field of frequency u> and amplitude 
A. Our first example is taken with generic values, A = 10, uj = 4, and the time run up to 
t = 4. In Fig. 21.1, we show how the system responds. The density barely changes from its 
ground-state value, but simply jiggles to left and right. The mean position begins to oscillate 
immediately, not with the frequency of the driving field, but rather with a period of about 
1/3. This corresponds to the lowest transition, since ei = 7r 2 /2 = 4.93, and = 4ei, so 
that W21 = 14.8, and the period is 2tt/uj ~ 0.4. Thus the time-dependent evolution of the 
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Figure 21.1: Electric field and expectation value of x for a particle in a box of length 1 subjected to the time-dependent field 
shown and given in the text. 

expectation value of operators contains information about transition frequencies. 

Next, we illustrate the concept of resonance. We run our calculation with a = 1, but 
sweep the driving frequency through u> = ujyi, as shown in Fig. 21.2. We see that, the nearer 
the driving frequency is to the transition frequency, the greater the effect of the driving field 
on the particle. Asw^ w 2 i, the perturbation grows indefinitely and pushes more and more 
energy into the particle. 

Another interesting regime of time-dependent problems is when the time variation is very 
slow. This is the adiabatic limit. For our particle in a box in an external electric field, if the 
external frequency u> <C w 12 , then the chances of exciting out of the ground-state in any small 
interval are exponentially small. Over long times, as the potential changes significantly from 
its original, these small chances add up, and the wavefunction does evolve. But if, instead 
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Figure 21.2: Expectation value of x and energy of a particle in a box of length 1 subjected to time-dependent fields of 
frequencies w — 13, 14, 15, and 16. 

of using the unperturbed eigenfunctions as our basis, we used the instantaneous eigenstates, 
i.e., those satisfying 

H{t) V n {t) = E n (t) (21.9) 

these chances would remain exponentially small, and so the system stays in the adiabatic 
ground state. Essentially, the potential is moving so slowly that the wavefunction of the 
particle simply rearranges itself to the new ground state at any given time. 

In Fig X, I've calculated the case where w = 1, t = 0.1, and a = 2000. Thus uj is much 
lower than the lowest transition, and uit <S 1. This yields a final electric field of a(wt) 2 /2 
or 10. The particle adiabatically follows the ground-state of the potential. This limit can be 
found by simply turning on a static external potential, and asking what the ground-state is. 

Exercise 88 Harmonic oscillator in a time-dependent electric field 

Consider a harmonic oscillator of frequency loq in an electric field. If it starts out at rest at 

the origin, show that its classical position is given by 

xM) = - f dt' sin(w (< - ?)) (21.10) 

Jo muj 

Then show that 

(f)(xt) = <p (x - x d (t)) exp(ice(t)) (21.11) 
is a solution to the time-dependent Schrddinger equation. How important is ct(t) ? 

In this framework, the time-dependent density is given by the usual formula of integrating 

the square of the wavefunction over all coordinates but one. The density operator may be 
written as 

N 

»(r) = E *(r - ii) (21.12) 

i=i 

and 

n (rt) = (*(a;)|n(r)|*(a;)) (21.13) 
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. For our particle in a box, it is simply 

CO 

n(xt) = \V(xt)\ 2 = £ c*(t) c k (t) ^(x) (j> k {x) (21.14) 

j,k=l 

A new quantity, that will play a central role in the development of TDDFT, is the current 
density. The operator is defined as 

j(r) = \ E (P^(r - n) + S(v - rj )Pi) (21.15) 

and its expectation value given by evaluation on the time-dependent wavefunction. For our 
particle in a box, 

j(xt) = Q(j)*(xt)d<t>(xt)/dx (21.16) 

Since, in any dx around a point x, the change in the number of electrons per unit time must 
equal the difference between the flow of electrons into dx from the left and out from the 
right, 

h(xt) = ~dj(xt)/dx (21.17) 

This is the equation of continuity. More generally, the time evolution of any (time-independent) 
operator is 

A = i[H(t),A] (21.18) 

Applying this to the density yields 

h(rt) = <tf(t)|i[f,n(r)]|tf(i)) = -Vj(i) (21.19) 

where the density commutes with all but the kinetic operator in the Hamiltonian, and its 
commutation (with a little algebra) yields the divergence of the current operator. 



21.2 Perturbation theory 

This concept is best seen in perturbation theory. Imagine integrating the equations of motion 
for the coefficients just a small amount in time, so that the excited states are only infinites- 
imally populated. Then one may insert the initial condition into the right-hand-side of Eq. 
X to find the probability of excitation. Next consider the effect of varying the external fre- 
quency, and one finds a huge response as it passes through a transition frequency. We find 
the transition rate, i.e., the probability of excitation per unit time, to be 

W(k) = P(k) = ... (21.20) 

which is known as Fermi's golden rule. 

The full range of behaviors of time-dependent systems is unnervingly large. In many cases, 
we are interested only in the response to a weak potential, i.e., 

H(t) = f + V + 5V(t) (21.21) 
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where 5V(t) is in some sense small. Usually we will consider potentials whose time-dependence 
is exp(iu!t + r/t) , where r\ — > 0, and the perturbation begins at — oo, and is fully turned on by 
t = 0. Under these conditions, we can ask what is the change in the system to first-order in 
the perurbation? More precisely, we want to know the change in observables of the system. 
Since this is a book on density functional theory, a key role will be played by the change in 
density in response to such a perturbation. We define the susceptibility of the system by 

Sp(xt) = J dx' J dt' x (x,x';t- f) 5v(x't') (21.22) 

i.e., \ tells you the density change throughout the system induced by any external potential. 
It is also known as the density-density response function and the non-local (meaning depends 
on x and x') susceptibility. Since we are doing DFT, we will write it more neatly as 

x[n Q ](x,x';t-t') =5p{xt)/5v{x't') (21.23) 

i.e., it is a functional derivative, and depends only on the ground-state density (for a non- 
degenerate system). Note that it is a function of t — t' alone, since the unperturbed system 
is static. Furthermore, we will speak of the retarded function, meaning x = for t' < t. We 
show in the rest of this section that x contains most of what we want to know about the 
response of a system. 

Next we specialize further to the specific case of the optical response of a system, i.e., to 
a long wavelength electric field. We want to find out the time-dependent dipole moment of 
our system, defined by 

d(t) = J dx p(xt)x (21.24) 

ie the first moment of the density. This is because the optical absorption is given by ... If we 

set v(x'ui) = Ex! , we find 

a{uj) = J dxxj dx' x'x(x, x'w) (21.25) 
for the dynamic polarizability, and 

a{u>) = 3a(wV/c (21.26) 

ie the moments of the susceptibility yield the complex dynamic susceptibility, which in turn 
give the optical response of the system. One can show this satisfies the well-known Thomas- 
Reich-Kuhn sum-rule, 

J™ckoa(uj) = 27 ^N (21.27) 

Thre is also the static sum-rule, which follows from the general Kramers-Kronig rule for the 
polarizability: 

2 f X — 9a(w) = 5Ra(0) (21.28) 

7T JO til 
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To see better what this means, we can write x m its Lehman representation, i.e., in terms 
of the levels of the system: 

X (xx'w) = MxWx') E - (21-29) 

We define: 

n jk (x) = (j\n(x)k) (21.30) 
the generalized transition dipole. Thus, taking the moments to construct a, we find 

x jk = J dx x n jk (x) (21.31) 
are the transition dipole moments for each transition. Then 

«M = E -P^ (21.32) 

To see the effect of specific transitions on the optical spectrum, we must use the well-known 
formula: 

— 1 — = V ( — 1 — ) + in 8(uj - loO) (21.33) 

UJ — UJ \ui — UJ J 

ie the real part is a principal value, and the imaginary part a delta-function. Thus 

a(u) = E 2cj j0 \x l0 \ 2 8(ij - lj q ) (21.34) 

The optical spectrum consists of a series of delta functions, whose intensity is determined by 
the dipole matrix elements. 

Exercise 89 Dipole matrix elements 

Calculate the dipole matrix elements for optical absorption for (a) the harmonic oscillator and 
(b) the particle in a box. 

We define the oscillator strength of a given transition as 

fji = 20^173 (21.35) 

These satisfy the TRK sum-rule, 

E fjjk = N (21.36) 

Note that all oscillator strengths from the ground-state are positive, but from an excited 
state, some are negative. 

Exercise 90 Oscillator strengths 

From above, now find oscillator strengths for both systems and show that they satisfy the 
Thomas-Reich-Kuhn sum-rule. 
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When discussing first-order perturbation theory, it is always useful to define response func- 
tions of the system at hand. A well-known one is the Green's function, which, for a one- 
electron system, may be written in operator form as 

G(z) = (z - H) { - 1] (21.37) 

where z is a complex number with positive imaginary part, and boundary conditions are 
chosen so that G(x,x') vanishes at large distances, just like the bound-states. In real space, 
G satisfies the differential equation: 

^V 2 + v(r) - z) G(rr'; z) = 5 (3) (r - br') (21.38) 



2 

Formally, one may write G in terms of sums over the eigenstates 

G(rr ' ;z) = E mj>^l. (21 . 39) 

i Z El 

Note that the Green's function is ill-defined whenever its argument matches that of an eigen- 
value. In such cases, we define 

G ± (rr', A) = G(rr', A ± irf) (21.40) 

where + refers to the retarded Green's function, while — refers to the advanced one. When 
Fourier-transformed, the Green's function yields the time-evolution of any initial state, for a 
time-independent Hamiltonian: 

tt(rt) = / dV f dt G+(rr', t - t') *(r'i') (21.41) 

since 

G+(t - = exp(-iH{t - t')) (21.42) 

Simple expression for G in terms of L and R phi. 

Relate density matrix to equal-time G 

Define \ and discuss meaning 

Relate % and G for non-interacting systems. 

Relate pair density to equal-time %■ 



21.3 Optical response 



Define polarizability, static and dynamic. 
Define oscillator strength. 
Define optical spectrum. 
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21.4 Questions about time-dependent quantum mechanics 

1. If a system is in a linear combination of two eigenstates, does its density change? 

2. Calculate the oscillator strengths for a particle in a box, and show that they satisfy the 
Thomas-Reich Kuhn sum rule. 

3. Calculate the imaginary part of the frequency-dependent susceptibility of the Id H-atom. 

4. Calculate the transmission and reflection amplitudes of the Id H-atom, assuming a par- 
ticle incident from the left. 

5. Calculate the retarded green's function for the Id H-atom. 

6. Using the first two levels of the particle in a box, calculate the Rabi frequency, and 
compare its time evolution with that of Fig X. 

7. Show how a harmonic oscillator evolves in a time-dependent electric field. 

8. Calculate the oscillator strengths for the harmonic oscillator, and show that they satisfy 
the Thomas-Reich Kuhn sum rule. 



Chapter 22 

Time-dependent density functional theory 



So far, we have restricted our attention to the ground state of nonrelativistic electrons in time- 
independent potentials. This is where most DFT research has been done, and has achieved its 
great successes, by providing a general tool for tackling electronic structure problems. In this 
chapter, we extend our view, to include time- dependent potentials. In subsequent chapters, 
we apply the methods of DFT to such potentials. 

22.1 Overview 

We will see below that, under certain quite general conditions, one can establish a one-to- 
one correspondence between time-dependent densities n(rt) and time-dependent one-body 
potentials v ex t(rt), for a given initial state. That is, a given evolution of the density can be 
generated by at most one time-dependent potential. This is the time-dependent analog of 
the Hohenberg-Kohn theorem, and is called the Runge-Gross theorem. Then one can define a 
fictitious system of non-interacting electrons moving in a Kohn-Sham potential, whose density 
is precisely that of the real system. The exchange-correlation potential (defined in the usual 
way) is then a functional of the entire history of the density, n(x), the initial interacting 
wavefunction ^(0), and the initial Kohn-Sham wavefunction, $(0). This functional is a very 
complex one, much more so than the ground-state case. Knowledge of it implies solution of 
all time-dependent Coulomb-interacting problems. 

An obvious and simple approximation is adiabatic LDA (ALDA), sometimes called time- 
dependent LDA, in which we use the ground-state potential of the uniform gas with that 
instantaneous and local density, i.e., v^ c [n](x) « w™ lf (n(x)). This gives us a working Kohn- 
Sham scheme, just as in the ground state. We can then apply this DFT technology to all 
problems of time-dependent electrons. These applications fall into three general categories: 
non-perturbative regimes, linear (and higher-order) response, and ground-state applications. 

The first of these involves atoms and molecules in intense laser fields, in which the field 
is so intense that perturbation theory does not apply. In these situations, the perturbing 
electric field is comparable to or much greater than the static electric field due to the nuclei. 
Experimental aims would be to enhance, e.g., the 27th harmonic, i.e., the response of the 
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system at 27 times the frequency of the perturbing electric field, or to cause a specific chemical 
reaction to occur (quantum control). Previously, only one and two electron systems could 
be handled computationally, as the full time-dependent wavefunction calculations are very 
demanding. Crude and unreliable approximations had to be made to tackle larger systems. 
But with the advent of TDDFT, which has become very popular in this community, larger 
systems with more electrons can now be tackled. 

When the perturbing field is weak, as in normal spectroscopic experiments, perturbation 
theory applies. Then, instead of needing knowledge of v xc for densities that are chang- 
ing significantly with time, we only need know this potential in the vicinity of the initial 
state, which we take to be a non-degenerate ground-state. These changes are character- 
ized by a new functional, the exchange-correlation kernel. While still a more complex beast 
than the ground-state exchange-correlation potential, the exchange-correlation kernel is much 
more manageable than the full time-dependent exchange-correlation potential, because it is a 
functional of the ground-state density alone. Analysis of the linear response then shows that 
it is (usually) dominated by the response of the ground-state Kohn-Sham system, but cor- 
rected by TDDFT via matrix elements of the exchange-correlation kernel. In the absence of 
Hartree-exchange-correlation effects, the allowed transitions are exactly those of the ground- 
state Kohn-Sham potential. But the presence of the kernel shifts the transition frequencies 
away from the Kohn-Sham values to the true values. Moreover, the intensities of the optical 
transitions can be extracted in the same calculations, so these also are affected by the kernel. 

Chemists and physicists have devised separate approaches to extracting excitations from 
TDDFT for atoms, molecules, and clusters. The chemists way is to very efficiently convert 
the search for poles of response functions into a large eigenvalue problem, in a space of 
the single-excitations of the system. The eigenvalues yield transition frequencies while the 
eigenvectors yield oscillator strengths. This allows use of many existing fast algorithms to 
extract the lowest few excitations. In this way, TDDFT has been programmed into most 
standard quantum chemical packages and, after a molecule's structure has been found, it is 
usually not too costly to extract its low-lying spectrum. Physicists, on the other hand, prefer 
to simply solve the time-dependent Kohn-Sham equations when a weak field has been turned 
on. Fourier-transform of the time-dependent dipole matrix element then yields the optical 
spectrum. Using either methodology, the number of these TDDFT response calculations 
for transition frequencies is growing exponentially at present. Overall, results tend to be 
fairly good (0.1 to 0.2 eV errors, typically), but little is understood about their reliability. 
Difficulties remain for the application of TDDFT to solids, because the present generation of 
approximate functionals (local and semi-local) lose important effects in the thermodynamic 
limit, but much work is currently in progress. 

The last class of application of TDDFT is, perhaps surprisingly, to the ground-state prob- 
lem. This is because one can extract the ground-state exchange-correlation energy from a 
response function, in the same fashion as perturbation theory yields expressions for ground- 
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state contributions in terms of sums over excited states. Thus any approximation for the 
exchange-correlation kernel yields an approximation to E xc . Such calculations are signifi- 
cantly more demanding than regular ground-state DFT calculations. Most importantly, this 
produces a natural method for incorporating time-dependent fluctuations in the exchange- 
correlation energy. In particular, as a system is pulled apart into fragments, this approach 
includes correlated fluctuations on the two separated pieces. While in principle all this is 
included in the exact ground-state functional, in practice TDDFT provides a natural method- 
ology for modeling these fluctuations. Simple ALDA approximations for polarizabilities of 
atoms and molecules allow Cg coefficients to be accurately calculated by this method, and 
attempts exist to generate entire molecular energy curves, from the covalent bond distance 
to dissociation, including van der Waals. Such calculations are much more demanding than 
simple self-consistent ground-state calculations, but may well be worth the payoff. 

22.2 Runge-Gross theorem 

The analog of the Hohenberg-Kohn theorem for time-dependent problems is the Runge-Gross 
theorem[?]. We consider N non-relativistic electrons, mutually interacting via the Coulomb 
repulsion, in a time-dependent external potential. The Runge-Gross theorem states that the 
densities n(x) and n'(x) evolving from a common initial state fy = ^(t = 0) under the 
influence of two potentials v ext (x) and v' ext (x) (both Taylor expandable about the initial time 
0) are always different provided that the potentials differ by more than a purely time-dependent 
(r-independent) function: 

A Vext {x) = v{x) - v'(x) ^ c (t) . (22.1) 

If so, there is a one-to-one mapping between densities and potentials, and one can construct 
a density functional theory. 

We prove this theorem by first showing that the corresponding current densities must differ. 
Because the Hamiltonians only differ in their external potentials, the equation of motion for 
the difference of the two current densities is: 

Aj(x)| t=0 = -i(* | [j(r), AH(t Q )] |*o) = -n (r)VA Uext (r, 0), (22.2) 

where ri (r) = n(r, 0) is the initial density. One can go further, by repeatedly using the 
equation of motion, and considering t = 0, to find [?] 

d k+1 Aj(x) = -no(r) Vd*Av ext (x) , (22.3) 

where dfi = (d k /dt k )\ t=0 , i.e., the k-th derivative evaluated at the initial time. If Eq. (22.1) 
holds, and the potentials are Taylor expandable about t , then there must be some finite k 
for which the right hand side of (22.2) does not vanish, so that 

Aj(.t) ^ 0. (22.4) 
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This is the first part of the Runge-Gross theorem, which establishes a one-to-one correspon- 
dence between current densities and external potentials. 

In the second part, we extend the proof to the densities. From continuity, 

An(i) = -V ■ Aj(x) (22.5) 

so that 

d k +2 An(x) = V ■ (n (r)Vd k Av ext (x)) . (22.6) 

Now, if not for the divergence on the right-hand-side, we would be done, i.e., if f(r) = 
dQAv ext (x) is non-zero for some k, then the density difference must be non-zero, since no(r) 
is non-zero everywhere. 

Does the divergence allow some escape from this conclusion? The answer is no, for any 
physical density. To see this, write 

/ d\ /(r) V ■ (n (r)V/(r)) = / d\ [V ■ (/(r)n (r) V/(r)) - n (r) | V/(r) | 2 ] (22.7) 

The first term on the right may be written as a surface integral at r = oc, and vanishes for 
any realistic potential (which fall off at least as fast as — 1/r). If we imagine that V/(r) 
is non-zero somewhere, then the second term on the right is definitely negative, so that the 
integral on the left cannot vanish, and its integrand must be non-zero somewhere. Thus 
there is no way for V/(r) to be non-zero, and yet have /V(n V/) vanish everywhere. 

Since the density determines the potential up to a time-dependent constant, the wavefunc- 
tion is in turn determined up to a time-dependent phase, which cancels out of the expectation 
value of any operator. 

22.3 Kohn-Sham equations 

We will see below that, under certain quite general conditions, one can establish a one-to- 
one correspondence between time-dependent densities n(rt) and time-dependent one-body 
potentials v ext (rt), for a given initial state. That is, a given history of the density can be 
generated by at most one history of the potential. This is the time-dependent analog of the 
Hohenberg-Kohn theorem. Then one can define a fictious system of non-interacting electrons 
that satisfy time-dependent Kohn-Sham equations 

i^-(rt) = + «s[n](x)j <j>j(x) . (22.8) 

whose density is 

n(x) = £ \^(x)\ 2 , (22.9) 
precisely that of the real system. We define the exchange-correlation potential via: 



v B {x) = v ext (x) + v H (x) + v xc (x) 



(22.10) 
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where the Hartree potential has the usual form, but for a time-dependent density. The 
exchange-correlation potential is then a functional of the entire history of the density, n(x), 
the initial interacting wavefunction ^(0), and the initial Kohn-Sham wavefunction, $(0). 
This functional is a very complex one, much more so than the ground-state case. Knowledge 
of it implies solution of all time-dependent Coulomb interacting problems. 

Exercise 91 Kohn-Sham potential for one electron 

By inverting the time-dependent Kohn-Sham equation, give an explicit expression for the 
time-dependent Kohn-Sham potential in terms of the density for one electron. 

In ground-state DFT, the exchange-correlation potential is the functional derivative of 
E xc [n]. It would be nice to find a functional of n(rt) for which v xc (rt) was the functional 
derivative, called the exchange-correlation action. However, as shown in Sec. X, this turns 
out not to be possible. 

22.4 Adiabatic approximation 

As we have noted, the exact exchange-correlation potential depends on the entire history of 
the density, as well as the initial wavefunctions of both the interacting and the Kohn-Sham 
systems: 

v xc [n; *(0), *(0)](rt) = « B [n;*(0)](rt) - v ext [n; *(0)](rt) - v B [n](rt). (22.11) 

However, in the special case of starting from a non-degenerate ground state (both interacting 
and non-), the initial wavefunctions themselves are functionals of the initial density, and so 
the initial-state dependence disappears. But the exchange-correlation potential at r and t 
has a functional dependence not just on n(r't) but on all n(v't') for < t' < t. We call this 
a history-dependence. 

The adiabatic approximation is one in which we ignore all dependence on the past, and 
allow only a dependence on the instantaneous density: 

<f[n](r f ) = <r°>W]M, (22.12) 

i.e., it approximates the functional as being local in time. If the time-dependent potential 
changes very slowly (adiabatically), this approximation will be valid. But the electrons will 
remain always in their instantaneous ground state. To make the adiabatic approximation 
exact for the only systems for which it can be exact, we require 

v^[n](rt) = < c M(r)U, l(r)=n(rt) (22.13) 

where vf£,[no](r) is the exact ground-state exchange-correlation potential of the density n (r). 
This is the precise analog of the argument made to determine the function used in LDA for 
the ground-state energy (see Sec. X). 
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In practice, the spatial locality of the functional is also approximated. As mentioned above, 
ALDA is the mother of all TDDFT approximations, but any ground-state functional, such 
as a GGA or hybrid, automatically yields an adiabatic approximation for use in TDDFT 
calculations. A simple way to go beyond the adiabatic approximation would be to include 
some dependence on, say, h(rt). 

Bureaucratically, ALDA should only work for systems with very small density gradients in 
space and time. It may seem like a drastic approximation to neglect all non-locality in time. 
However, we've already seen how well LDA works beyond its obvious range of validity for the 
ground-state problem, due to satisfaction of sum rules, etc., and how its success can be built 
upon with GGA's and hybrids. We'll soon see if some of the same magic applies to TDDFT. 

Exercise 92 Continuity in ALDA 

Does a Kohn-Sham ALDA calculation satisfy continuity? Prove your answer. 

22.5 Questions on general principles of TDDFT 

1. If a system is in a linear combination of two eigenstates, does its density change? 

2. In a time-dependent Kohn-Sham calculation starting from the ground state, will electronic 
transitions occur at frequencies that are differences between ground-state Kohn-Sham 
orbital energies? 

3. Is there anything special about the response of a particle in a harmonic well to an external 
electic field? 

4. Does the Runge-Gross theorem guarantee that the current density in a TDDFT calcula- 
tion is correct? 

5. Can more than one potential produce the same time-dependent density? 

6. If a system is subjected to a driving potential that, after a time, returns its density to its 
initial value, will the potential return to its initial value? 

7. Repeat question above, but replace density with wavefunction. 

8. Repeat both questions above, within the adiabatic approximation. What happens if your 
functional has some memory? 

9. Consider a He atom sitting in its ground-state. Imagine you have a time-dependent 
potential that pushes it into its first excited state (ls2,s). Now imagine the exact Kohn- 
Sham calculation. Describe the final situation. 
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Chapter 23 
Linear response 



Note that the difference between n(x) and n'{x) is non-vanishing already in first order of 
v(x) — v'(x), ensuring the invertibility of the linear response operators of section ??. 

23.1 Dyson-like response equation and the kernel 

When the perturbing field is weak, as in normal spectroscopic experiments, perturbation 
theory applies. Then, instead of needing knowledge of i> xc for densities that are changing 
significantly with time, we only need know this potential in the vicinity of the initial state, 
which we take to be a non-degenerate ground-state. Writing n(x) = n(r) + 5n(x), we have 

v xc [n + 5n](x) = v xc [n](r) + / dt' f d 3 r' / x0 [ra](rr', t - t') 5n(r't'), (23.1) 
where / xc is called the exchange-correlation kernel, evaluated on the ground-state density: 

f xc [n](rr',t- t') = 5v xc (rt)/5n(r'tf). (23.2) 

While still a more complex beast than the ground-state exchange-correlation potential, the 
exchange-correlation kernel is much more manageable than the full time-dependent exchange- 
correlation potential, because it is a functional of the ground-state density alone. To un- 
derstand why f xc is important for linear response, we define the point-wise susceptibility 
x[no](rv', t — t') as the response of the ground state to a small change in the external poten- 
tial: 

5n(x) = J dt' J d 3 r' X [n }(rr',t- t') 8v ext {r't'), (23.3) 

i.e., if you make a small change in the external potential at point r' and time t', \ tells you 
how the density will change at point r and later time t. Now, the ground-state Kohn-Sham 
system has its own analog of \, which we denote by \s, which says how the non-interacting 
KS electrons would respond to Sv s (r'f), which is quite different from the interacting case. 
But both must yield the same density response, so using the definition of the KS potential, 
and of the exchange-correlation kernel, we find the key equation of TDDFT linear response: 

X(rr'w) = Xs(rr'w) + / cPn f d 3 r 2 Xs(rriw) ( j^r^j + /xc(rir 2 o;)j x(r 2 r'w), (23.4) 
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where all objects are functional of the ground-state density. This is called a Dyson-like 
equation, because it has the same mathematical form as the Dyson equation relating the 
one-particle Green's function to its free counterpart, with a kernel that is the one-body 
potential. However, it is an oddity of DFT that it has this form, since x is a two-body 
response function. 

This equation contains the key to electronic excitations via TDDFT. This is because when 
uj matches a true transition frequency of the system, the response function % blows up, i.e., 
has a pole as a function of uj. Thus x s has a set of such poles, at the single-particle excitations 
of the KS system: 

*(-^) = 2 £ ^y^^') .^^^ (23 . 5) 

ioce,aunocc ^ V-i ^-a) ' ^+ 

In the absence of Hartree-exchange-correlation effects, \ = Xs, and so the allowed transi- 
tions are exactly those of the ground-state KS potential. But the presence of the kernel in 
Eq. (23.4) shifts the transitions away from the KS values to the true values. Moreover, 
the strengths of the poles can be simply related to optical absorption intensities (oscillator 
strengths) and so these also are affected by the kernel. 

In linear response, this implies that the exchange-correlation kernel has the form: 

fit (rr', t-t')= Sj f^f 1 U)=n(r t) S(t - t>). (23.6) 

When Fourier-transformed, this implies /^ la is frequency-independent. 



23.2 Casida's equations 

The various methods for extracting excitations have been described in the overview. Here we 
focus on the results. Casida showed that, finding the poles of \ is equivalent to solving the 
eigenvalue problem: 

J2Q qq ,(uj)v q , = Qv q , (23.7) 
«' 

where q is a double index, representing a transition from occupied KS orbital i to unoccupied 
KS orbital a, uj q = e a — e u Q = uj 2 , and ( & 5 (r) = 0*(r)0 a (r). The matrix is 

Q qq ,(oj) = 5 qq ,Q. q + 2^7 q {q\f HXC (uj)\q'), (23.8) 

where 

( g |/ H xcHk'> = fd\ /d'V$;(r)/ EXC (r,r»$,,(r'). (23.9) 

In this equation, / HXC is the Hartree-exchange-correlation kernel, l/|r — r'| + / xc (r, r',a;). 
A selfconsistent solution of Eq. (23.7) yields the excitation energies co and the oscillator 
strengths can be obtained from the eigenvectors [?]. In the special case of an adiabatic 
approximation, these equations are a straightforward matrix equation. Algorithms exist for 
extracting just the lowest transitions. 



23.3. SINGLE-POLE APPROXIMATION 

23.3 Single-pole approximation 
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Chapter 24 
Performance 



24.1 Sources of error 

24.2 Poor potentials 

A problem was noticed early on. Our favorite approximations to the exchange-correlation 
energy yield fairly lousy exchange-correlation potentials. These potentials are especially un- 
pleasant at large distances from Coulombic systems. They are too shallow, making Koopman's 
theorem fail drastically. This means that all the higher states are too shallow. In fact the 
Rydberg series, i.e., an infinite sequence of excitations as E — > from below, characteristic 
of the — 1/r decay, does not exist at all. How can one correct KS transition frequencies if 
they do not exist? 

The standard answer has been to asymptotically correct the potentials, to make them mimic 
the true KS potential more accurately. A variety of schemes exist that do this to varying 
degrees of accuracy. However, we note that all such schemes then produce a potential that 
is not a functional derivative of an exchange-correlation energy. 

24.3 Transition frequencies 

In any event, for many systems of real interest, it is the low-lying valence excitations that are 
of interest, and in these regions, the potentials are quite accurate. Table X shows... 

Hybrid functionals in TDDFT: An important point to note here concerns imple- 
mentation of hybrid functionals in TDDFT. There are two ways a hybrid might be coded, 
which in the case of excitations, can yield very different results. In the first, 'correct' DFT 
approach, the optimized effective potential method should be used, to guarantee that the 
effective potential is always a local, multiplicative one, as required by Kohn-Sham theory. 
However, in practice, Hartree-Fock calculations are much easier, and have been around since 
time immemorial. Thus the potential produced in many quantum chemical codes consists 
of a fraction of the HF non-local potential. This has negligible effect for ground-state and 
occupied orbitals, but a large effect on the unoccupied levels. 
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At this point, we have surveyed all that one needs to perform a TDDFT calculation 
of excitation energies using a modern code, either a quatum chemistry one or a real-time 
evolution. The remainder of this chapter, and all of the next, are devoted to delving deeper 
into understanding how TDDFT works and what are its current limitations. 

In this section, we report calculations for the He and Be atoms using the exact ground- 
state Kohn-Sham potentials, the three approximations to the kernel mentioned in the previous 
section, and including many bound-state poles in Eq. (??), but neglecting the continuum. 
The technical details are given in Ref. [?]. Table I lists the results, which are compared with 
a highly accurate nonrelativistic variational calculations[?, ?] In each symmetry class (s, p, 
and d), up to 38 virtual states were calculated. The errors reported are absolute deviations 
from the exact values. The second error under the Be atom excludes the 2s — > 2p transition, 
for reasons discussed in the next section. 

The effect of neglecting continuum states in these calculations has been investigated by van 
Gisbergen et al.[?], who performed ALDA calculations from the exact Kohn-Sham potential in 
a localized basis set. These calculations were done including first only bound states, yielding 
results identical to those presented here, and then including all positive energy orbitals allowed 
by their basis set. They found significant improvement in He singlet-singlet excitations, 
especially for Is — > 2s and Is — > 3s. Other excitations barely changed. Assuming inclusion 
of the continuum affects results with other approximate kernels similarly, these results do not 
change the basic reasoning and conclusions presented below, but suggest that calculations 
including the continuum may prove to be more accurate than those presented here. 

24.4 Atoms 

24.5 Molecules 



24.6 Strong fields 
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Table 24.1: Singlet/triplet excitation energies in the helium and beryllium atoms, calculated from the exact Kohn-Sham 
potential by using approximate xc kernels (in millihar trees), and using the lowest 34 unoccupied orbitals of s and p symmetry 
for He, and the lowest 38 unoccupied orbitals of s, p, and d symmetry for Be. Exact values from Ref. [?] for He and from 
Ref. [?] for Be. 



Singlet/triplet shifts 





UKS 


ALDA 


X 


SIC 


hybrid 


exact 


Transitions from the Is state in He atom 


2s 


746.0 


22/-11 


20/-25 


19/-16 


14/-19 


12/-18 


3s 


839.2 


6.9/-2.4 


5.8/-4.9 


5.6/-3.6 


5.1/-4.6 


3.3/-4.2 


4s 


868.8 


3.1/-0.9 


2.5/-1.7 


2.4/-1.3 


2.4/-2.0 


1.3/-1.6 


5s 


881.9 


1.6/-0.4 


1.3/-0.8 


1.3/-0.6 


1.3/-1.0 


0.6/-0.8 


6s 


888.8 


1.0/-0.3 


0.8/-0.5 


0.7/-0.4 


0.8/-0.7 


0.4/-0.5 


2p 


777.2 


-0.8/-7.4 


7.2/-8.4 


6.1/ 0.2 


2J/-3.4 


2.7/-6.6 


3p 


847.6 


0.7/-1.9 


2.5/-2.3 


2.2/-0.5 


1.4/-1.3 


1.0/-2.0 


Ap 


872.2 


0.4/-0.7 


1.1/-0.9 


1.0/-0.2 


0.6/-0.5 


0.5/-0.8 


5p 


883.6 


0.2/-0.4 


0.6/-0.5 


0.5/-0.1 


0.4/-0.3 


0.2/-0.4 


6p 


889.8 


0.1/-0.3 


0.3/-0.3 


0.3/-0.1 


0.2/-0.2 


0.1/-0.3 


err 


57 


32 


31 


29 


12 




Transitions from the 2s state in Be atom 


3s 


244.4 


7.1/-5.7 


10.9/-10.6 


10.3/-1.4 


6.6/-4.6 


4.7/-7.1 


4s 


295.9 


2.5/-1.6 


3.6/-2.5 


3.5/-0.6 


2.6/-1.6 


1.4/-2.0 


5s 


315.3 


1.1/-0.7 


1.7/-1.0 


1.6/-0.3 


1.2/-0.7 


0.6/-0.9 


6s 


324.7 


0.6/-0.4 


0.9/-0.5 


0.9/-0.2 


0.7/-0.4 


0.3/-0.5 


2p 


132.7 


56/-42 


55/-133 


53/-53 


10/-88 


61/-32 


3p 


269.4 


2.0/-4.3 


6.4/-4.2 


5.6/ 1.1 


4.2/-2.8 


4.8/-1.5 


Ap 


304.6 


0.3/-1.4 


2.1/-1.2 


1.9/ 0.3 


1.3/-0.8 


1.7/-4.1 


5p 


319.3 


0.1/-0.6 


1.0/-0.5 


0.9/ 0.1 


0.6/-0.2 


0.2/ 0.0 


6p 


326.9 


0.0/-0.4 


0.5/-0.3 


0.4/ 0.0 


0.3/-0.2 


0.1/-0.1 


3d 


283.3 


-5.4/-2.8 


1.8/-2.0 


0.9/ 3.2 


-1.4/ 1.2 


10.3/-0.6 


id 


309.8 


-1.4/-1.1 


0.8/-0.9 


0.5/ 0.9 


-0.2/ 0.6 


3.6/-0.2 


err 


138 


56 


144 


73 


136 




err' 


45 


41 


37 


44 


29 
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25.1 Currents 

25.2 Initial-state dependence 

Not much attention is paid in the literature to the initial-state dependence in the Runge-Gross 
theorem, and its often not mentioned explicitly. But TDDFT would be useless if we always 
had to account for it, since we'd have a different functional for every initial wavefunction. 

In practice, this is not really a problem. Often, one restricts oneself to starting in a 
non-degenerate ground state. In that case, by virtue of the ground-state Hohenberg-Kohn 
theorem, the initial wavefunction is an (implicit) density functional, and so the entire potential 
is a density functional. 

In fact, one can be even more general, and include any initial wavefunction that can be 
generated by some time-evolution of the system from a non-degenerate ground-state. Suppose 
one can find such a pre-history, using an external potential v(x), with t running from — tp to 
0. In that case, simply tack on this pre-history, and apply TDDFT with the initial ground- 
state functional, begining from — tp. The original initial potential (at t = 0) differs from 
what it would be if we began in the ground-state, because of the memory-dependence of the 
functional. 

Can all initial-states be generated by a pseduo-prehistory begining in a ground-state? This 
seems unlikely, since one only has the freedom to vary v ext (rt), but the wavefunction is a 
function of many variables. But presumably most of physical interest are of this kind. 

We will not discuss further the initial-state dependence, and assume from here on that the 
initial state is a non-degenerate ground state. 

25.3 Lights, camera, and... Action 

Before discussing the action in time-dependent DFT, I first point out that, for TDDFT, there 
is no analog of the ground-state energy. That is, there is no one functional in which we have 
an overarching interest. In a time-dependent problem, we wish to know the time-dependent 
expectation value of the density, which is determined by the Kohn-Sham potential. But in 
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the ground-state problem, what we really care about is the energy, which will determine the 
geometry of our molecules and solids. 

Nevertheless, the action is the mathematically analogous quantity, defined as 

A[%= J^Ut {^W -i^m)) (25.1) 
The principle of least action says that 

is satisfied by the wavefunction obeying the time-dependent Schrodinger equation. It would 
be very nice to define an exact A[n] whose variations with respect to the density yielded the 
time-dependent Schroedinger equation. 

The original Runge-Gross paper defined this action simply as 

A[n] = A[f[n]} (25.3) 

However, this has since been shown to lead to several inconsistencies. Essentially, variations 
in the density restrict variations in the wavefunction to only those that are t'-representable, 
i.e., form 



25.4 Solids 

25.5 Back to the ground state 

The last class of application of TDDFT is, perhaps surprisingly, to the ground-state problem. 
This is because one can extract the ground-state exchange-correlation energy from a response 
function, in the same fashion as perturbation theory yields expressions for ground-state con- 
tributions in terms of sums over excited states. To do this, we simply note that the equal-time 
susceptibility, xi^'^ — f = 0), yields the density-density correlation function, i.e., the pair 
density of chapter X. Thus 

r4(r, r') = - X \rr', t -> f )/n(r) - ^ 3 >(r - r'), (25.4) 

where we have now generalized Eq. (23.4) to arbitrary A by using a kernel A/|ri — r 2 | + 
/x C (rir 2 cj). Inserting this expression for the hole into the adiabatic connection formula for 
the energy, we find: 

^-/cfr/^jff .... (25.5) 

giving the exchange-correlation energy directly in terms of \ X - Any approximation for / xc , 
whose A-dependence can be simply extracted via scaling, yields an approximation to S xc . 

Most importantly, Eq. (25.5) produces a natural method for incorporating time-dependent 
fluctuations in the exchange-correlation energy. In particular, as a system is pulled apart 
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into fragments, Eq. (25.5) includes correlated fluctuations on the two separated pieces. 
While in principle all this is included in the exact functional, in practice TDDFT provides 
a natural methodology for modelling these fluctuations. Simple ALDA approximations for 
polarizabilities of atoms and molecules allow Cq coefficients to be accurately calculated by 
this method, and attempts exist to generate entire molecular enery curves, from the covalent 
bond distance to dissociation, including van der Waals. Such calculations are much more 
demanding than simple self-consistent ground-state calculations, but may well be worth the 
payoff. 

25.6 Multiple excitations 



25.7 Exact conditions 
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Appendix A 

Math background 



A.l Lagrange multipliers 

To illustrate how Lagrange multipliers work, we do a simple example on a function of two 
variables, instead of functionals. 

Suppose we wish to maximize the function f(x,y) = x + y subject to the constraint 
x 2 + y 2 = 1. The inelegant approach is to simply solve the constraint equation for y as a 
function of x, finding y = \/l — x 2 , and substitute into /, yielding 

f(x) = x + Vl-x 2 (A.l) 

Then we find the extrema of / by differentiation: 

f = 1 - (A.2) 
dx y/1 - x 2 ' 

which we set to zero, finding x = l/v2. Then y = 1/V2 also, yielding the maximum 

f = V2. 

But a much more elegant way is to introduce a Lagrange multiplier. Write the constraint 
as g(x, y) = 0, where g = x 2 + y 2 — 1. Then we define 

h = f-ng (A.3) 

where fj, is a constant, independent of x and y. We then minimize h freely, with no restriction 
between x and y. Write 

dh dh „ ,. ,. 

- = l-2,.x g- y =l-2,y (A.4) 

and set both to zero, yielding (x,y) = (l,l)/(2/i). Then we enforce the condition g = 0, 
yielding fj, = ±l/v2, and choose the positive sign for a maximum. This method avoids 
differentiating square roots, etc. This is especially important if we cannot solve the constraint 
equation analytically. 
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A.2 Properties of the 5-function 

We review briefly the useful properties of 5-f unction. The (^-function is an infinitely narrow, 
infinitely tall peak at x = 0, chosen so that its area is exactly 1. Thus we can write several 
limiting formulas: 

6{x) = lim (6(.x + a/2) - Q(x - a/2)) /a, 
= lim — = exp (—X 2 hi 2 ) 

= lim- 9 7 - (A.5) 

where 6(x) is the step function (= for x < 0, 1 for x > 0), and the 5-function looks like 
(with a finite width) what's shown in Fig X. 

The usefulness of the S function is that, for any reasonably smooth function: 

Jdxf(x)S(x) = f(Q) (A.6) 

assuming the integration interval includes the origin. This can be used to extract the value 
of the function at any point x Q : 

Jdxf(x)5(x-x ) = f(x ) (A.7) 

again assuming the integration interval includes the point x . 
Another handy formula is the integral over the 5-function: 

T dx' 5(x') = 9(x) (A.8) 

which is just the step function, i.e., you get zero unless the interval includes the (^-function, 
when you get 1. Another way to see this is that the first formula of Eq. (A.5) implies that 
the S function is the derivative of the step function. 
The Fourier transform of the J-function is 

dp 



S(x) = I — - exp(ipx) (A. 9) 



i.e., it's just a constant. 
Finally, one can show: 



where f'(x) = df /dx. 
Lastly, with some care, 



'(/(*))= £ %uf (A.10) 

i,/(x,)=0 1/ v x >)\ 



1 5'(x)f(x) = -/'(0) (A.ll) 

assuming the interval includes 0. 

In fact, its not strictly a function at all, but properly called a distribution. 
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A. 3 Fourier transforms 
Include Kramers-Kronig relation. 



Appendix B 

Results for simple one-electron systems 



B.l Id H atom 

The Hamiltonian for this problem is 

Integrating the Schrodinger equation through x = yields the cusp condition: 

0'ft = -Z<p(Q) (B.2) 

The ground-state energy is just e = —Z 2 /2, and T = —e , V = 2e Q . The normalized orbital 
is: 

exp(-Z|x|) (B.3) 

and ground-state density 

n(x) = Z exp(-2Z\x\) (B.4) 

We can also extract continuum excited states, which are useful for TDDFT. With the usual 
scattering boundary conditions: 

4>s{x) = exp(ifcx) + r exp(— ikx), x<0 

= texp(ikx), x>Q (B-5) 

where r and t are given by: 

iZ ik „. 

r = ~, Ti, t=- r B.6) 

k-iZ' Z + ik v ' 



and k = v2e, E > 0. Furthermore, when H is given by (B.l), the solution to: 

(e + -H)g > (x,x';e) = 5{x - x') (B.7) 
is the outgoing Green's function: 

g>(x, x': e) = e^~ x I - ^= . (B.; 

y V ' 7 iv^ ^ iV2e + Z ' 
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and in terms of g > , the susceptibility is 

X{x, x'\ e) = \jn{x)n{x') [g > (x, 2'; e + e ) + ff > *(a;, 2', -e + e )] (B.9) 
From this we can calculate both the static and dynamic polarizabilites: 

«(0) = - j|j (B.10) 

M«M) = (B.11) 

where Q = w/\eo\ 



B.2. HARMONIC OSCILLATOR 

B.2 Harmonic oscillator 
For the Hamiltonian 

the energies are 



H 



1<P 1 2 2 

-2d^ + 2 UJX 



n + -jw, n = 0,1,2,..., 
The lengthscale is given by x Q = y^, and the wavefunctions are 

1 



4>n{x) 



x ^2 n n\ 



\x 



(B.12) 
(B.13) 

(B.14) 



where H n (y) are the Hermite polynomials. The ground state and first 2 excited states wave 
functions are 

1 

_.x 0v /7r_ 



2 

.x 0v /7r 



-if.y 

1 



■I'D 



i 



(B.15) 
(B.16) 
(B.17) 
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B.3 H atom 
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Appendix C 
Green's functions 
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Appendix D 
Further reading 



Chapter 1 

There are now a variety of sources to choose from to learn about density functional theory 
in general. 

• Perhaps the best overview comes from a recent summer school, A Primer in Density 
Functional Theory, ed. C. Fiolhais, F. Nogueira, and M. Marques (Springer-Verlag, NY, 
2003). 

• A good general-purpose book for solid-state physicists about electronic structure is 
Richard Martin's book. Another is Vignale's. 

• For chemists, a very useful guide is A chemist's guide to density functional theory, W. 
Koch and M.C. Holthausen, (Wiley-VCH, Weinheim, 2000). Part A deals with theory, 
while part B has a very useful survey of properties and how well approximations do for 
them. 

• A standard in physics for many years has been: R.M. Dreizler and E.K.U. Gross, Density 
Functional Theory (Springer-Verlag, Berlin, 1990). This is very careful and rigorous, 
but uses some higher-level physics concepts. 

• In chemistry, we have Density Functional Theory of Atoms and Molecules, R.G. Parr 
and W. Yang (Oxford, New York, 1989), which includes many concepts particularly useful 
for chemistry, such as reactivitity theory. 

• A nice discussion of using Kohn-Sham orbitals to understand chemistry appears in A 

quantum chemical view of density functional theory, E.J. Baerends and O.V. Gritsenko, 
J. Phys. Chem. A 101, 5383 (1997). 

Chapter 3 

There are many books on basic quantum mechanics. 

• The one I enjoy the most is Griffiths.. 
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A great old one that demonstrates the unity of quantum is Morisson... 
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Discussion of questions 



Chapter 1 

1. Why is a Kohn-Sham calculation much faster than a traditional wavefunction calcu- 
lation? 

Because a Kohn-Sham calculation requires solving only a one-electron problem (self- 
consistently) and occupying N levels, rather than solving a coupled differential equation 
for a function of 3N variables. 

2. If you evaluate the kinetic energy of a Kohn-Sham system, is it equal to the physical 
kinetic energy? 

No. The kinetic energy operator is the same for the two systems, but their wavefunctions 
are different. 

3. Repeat above question for {l/r), which can be measured in scattering experiments. 

Yes. (1/r) = / d 3 r n(v)/r, and so depends only on the one-electron density, which is 
defined to be the same in both systems. 

4. Why is the density cubed in the local approximation for T s ? (see section ?? for the 
answer). 

Performing a dimensional analysis, we note that n ~ 1/L in one dimension, and since 

T ~ 1/L 2 = Jdxn IJ (x), p must be 3. 

5. Suppose you'd been told that T}° c is proportional to J dxn 3 (x), but were not told the 
constant of proportionality. If you can do any calculation you like, what procedure 
might you use to determine the constant in T s loc ? 

A physicist might answer that, since T s becomes more and more accurate as the number 
of electrons in the box grows, take the limit asiV^oc. We saw in the chapter how that 
yields the 7r 2 /6 constant. On the other hand, a practical chemist might answer that, if 
all I want to do is solve problems that look like one electron in a box, I should simply fit 
the constant to the value for that problem. 
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Both answers have value, but reflect different priorities. The chemist's answer is likely 
to be more accurate than the physicist's for problems that are similar to the fitted one. 
Errors are likely to be both positive or negative. On the other hand, the physicist's answer 
is likely to be more accurate over a broad range of problems, with systematic errors (but 
larger than the chemist's on the fitted systems). 

6. In what way will T s loc change if spin is included, e.g., for two electrons of opposite 
spin in a box? 

A detailed answer is given in section 9.2. Here we just doubly occupy each level in the 
box. Thus. ..(Rudy, fill in). 

7. // we add an infinitesimal to n(x) at a point, i.e., eS(x) as e — > 0, how does T g loc 
change? 

Again, a detailed answer is given in section 2.2. But we can figure it out easily in this 
case. Rudy... 

8. The simplest density functional approximation is a local one. What form might a cor- 
rection to Tl° c take? 

A detailed answer is given in section 17.2. But we can make intelligent suggestions here. 
The most obvious is that, since using the density at a point works so robustly, surely 
including information on the gradient of the density can help accuracy. Rudy... 

9. In the chapter and exercises, Tj| oc does very well. But there are cases where it fails 
quite badly. Can you find one, and say why? 

Consider a particle in a large box, at whose center is a tall but finite barrier. Ignoring the 
exponentially small tunneling contribution, the lowest molecular orbital is 

<f>(x) = (Mx) + Mx))/V2 (E.l) 

where 4>a(x) is the ground-state wavefunction for the electron to be in the left side, 
and 4>b{x) is that for the right. The exact kinetic energy is then \(Ta + T B ) = Ta by 
symmetry. On the other hand, our local functional uses n(x) = {tia{x) + n B (x))/2, so 
that T s loc = (ri oc + T|j oc )/8 = 7^/4, i.e., it is about 4 times too small! 

10. Does T s loc [n*] work for a single electron in an excited state, with density n*(r)? 

For ANY single excited state of the particle in the box, T s loc [n*] = T s loc [n], so it fails very 
badly in this form. Need to show explicit failure for particle in box. Rudy... 

However, Rudy points out that, since our electrons are not interacting, he can find the 
energy of the first excited state by writing n*{x) = n 2 {x) — n\{x), where ni(x) is the 
density with i particle in the box. Then T[n*} = T\ oc - Tj 00 = 17.7/L 2 



205 



Chapter 2 

1. Compare the functional derivative of T g vw [n] with T s ! oc [n] for some sample one- 
electron problem. Comment. 

You might take either the particle in the box or the harmonic oscillator or the one- 
dimensional hydrogen atom. You should find that a plot of the two functional derivatives 
gives you two very different curves in each case. To understand this, think of functions. 
Imagine one function being a horizontal line, and a second function having weak but 
rapid oscillations around that line. The second function is a good approximation to first 
everywhere, but its derivative is very different. 

2. Devise a method for deducing if a functional is local or not. 

Take the second functional derivative. If its proportional to a delta function, the function 
is local. Note that just saying if the integrand depends only on the argument at r is not 
enough, since integrands can change within functional, e.g., by integrating by parts. 

3. Is there a simple relationship between T s and / dx n(x)5T s /Sn(x)? 

You should have found that the integral is three times the functional for the local ap- 
proximation, but equal to the functional for von Weisacker. 

4. For fixed particle number, is there any indeterminancy in the functional derivative 
of a density functional? 

If the particle number does not change, then / d 3 r Sn(r) = 0. Thus addition of any con- 
stant to a functional derivative does not alter the result, so that the functional derivative 
is not determined up to a constant. 
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Chapter 3 

1. Suggest a good trial wavefunction for a potential that consists of a negative delta 
function in the middle of a box of width L. 

The simplest guess is just that of the particle in the box, ^2/Lcos(nx/L). Its energy is 
obviously 7r 2 /2L 2 — 2Z/L, where Z is the strength of the delta function. Much better 
(in fact, exact) is to take linear combinations of the delta-function solutions, exp(— a\x\) 
and exp(a|x|) vanishing at the edges. 

2. What is the effect of having nuclear charge Z =f 1 for the 1-d H-atom? 

It alters the length scale of the wavefunction without changing its normalization, i.e., 
\/~Z exp(— Z\x\). We see shortly how this can be a handy trick. 

3. What is the exact kinetic energy density functional for one electron in one- dimension? 
It is the von Weisacker functional, 

T * VW = \ II dx l^' 2 = \ II dx l"'M| 2 Mz)- (E.2) 

4. Is the HF estimate of the ionization potential for 1-d He an overestimate or under- 
estimate? 

The ionization potential is / = E\ — E 2 for a two-electron system, where E N is the 
energy with N electrons. Since _E HF > E, with equality for one electron, 7 HF < I. 

5. Consider the approximate HF calculation given in section 4.2. Comment on what it 
does right, and what it does wrong. Suggest a simple improvement. 

The HF calculation correctly writes the wavefunction as a spatially symmetric singlet, but 
incorrectly approximates that as a product of separate functions of x\ and x 2 (orbitals). 
It correctly has satisfies the cusp condition at the nucleus, and changes decay at large 
distance, unlike the scaled orbital solution. But its energy is not low enough. 

6. Which is bigger, the kinetic energy of the true wavefunction or that of the HF wave- 
function for Id He? Hooke's atom? 

The virial theorem applies to this problem (see section X), and, since all potentials are 
homogeneous of order— 1, E = —T = —V/2, where V includes both the external and 
the electron-electron interaction. Since _E HF > E, then T HF < T. For Hooke's atom, 
the situation is more complicated, since there are two different powers in the potentials. 

Chapters 4-6 

1. The ground-state wavefunction is that normalized, antisymmetric wavefunction that has 
density n(v) and minimizes the kinetic plus Coulomb repulsion operators. 
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2. The Kohn-Sham wavefu notion of density n(r) is that normalized, antisymmetric wave- 
function that has density n(r) and minimizes the kinetic operator. 

3. The Kohn-Sham kinetic energy is not |/ dx |0'j| 2 + ^l 2 , because these orbitals could 
have come from anywhere. There is no reason to think they are Kohn-Sham orbitals. 
They might be, e.g., Hartree-Fock orbitals. The way to get T s is to construct the density 
n{x) and then find v s (x), a local potential that, with two occupied orbitals, yields that 
density. The kinetic energy of those two orbitals is then T s . Again, if we alter one 
of the original orbitals, the change in T s is not the kinetic contribution directly due to 
that orbital. Rather, one must calculate the new density, construct the new Kohn-Sham 
potential, get the kinetic energy of its orbitals, and that tells you the change. 

All the same reasoning applies for -E x [?i]. 

4. See above. One can do a HF calculation, yielding orbitals and density. Then E^ F is 
simply the Fock integral of those orbitals. But this is not _E x [n HF ], which could only be 
found by finding the local potential v s (r) whose orbitals add up to n HF (r), and evaluating 
the Fock integral on its orbitals. 

5. The additional flexibility of spin-DFT over DFT means that its much easier to make good 
approximations for spin-polarized systems. 

6. Given v S!J (r), solve the Kohn-Sham equations for up and down non-interacting electrons 
in these potentials, and construct the total density. Then ask what single potential all 
electrons must feel in the Kohn-Sham equations to reproduce that density. Obviously, 
they coincide when the system is unpolarized, so that v sa (r) = v s (r), and when fully 
polarized, i> S f(r) = v s (r), v S j(r) = or undetermined. 

7. No. See further work whead. 



